A population model for the 2017/18 listeriosis outbreak in South Africa

We introduce a compartmental model of ordinary differential equations for the population dynamics of listeriosis, and we derive a model for analysing a listeriosis outbreak. The model explicitly accommodates neonatal infections. Similarly as is common in cholera modeling, we include a compartment to represent the reservoir of bacteria. We also include a compartment to represent the incubation phase. For the 2017/18 listeriosis outbreak that happened in South Africa, we calculate the time pattern and intensity of the force of infection, and we determine numerical values for some of the parameters in the model. The model is calibrated using South African data, together with existing data in the open literature not necessarily from South Africa. We make projections on the future outlook of the epidemiology of the disease and the possibility of eradication.


Introduction
The World Health Organisation (WHO, 2018) reported [1] the 2017/2018 outbreak of listeriosis in South Africa as the largest global outbreak recorded to date. Between 1 January 2017 and 14 March 2018, the National Institute for Communicable Diseases (NICD) reported 978 confirmed cases of whom 42% were neonates, infected via vertical transmission from infected pregnant mothers. Listeriosis is transmitted through the contamination of infected food products with a Gram-positive non-spore forming, facultatively anaerobic bacillus called Listeria monocytogenes [2] (Lamont et al 2011). The outbreak was predominantly attributed to sequence type 6 (ST6) while other sequence types attributed to 9% of infections [1] (WHO, 2018). Listeriosis is reported to occur 18 times more frequently in pregnant than non-pregnant women [3] (Jackson et al., 2010). Infection is usually asymptomatic or may primarily manifest as a non-specific flu-like illness which is often overlooked, particularly during pregnancy, leaving the mother and her unborn infant at risk of serious illness. The infection may be vertically transmitted while the baby is in utero, either by translocation across the placenta or by aspiration of contaminated amniotic fluid or or during delivery of the infant. Foetal infection (during the 1st and 2nd trimester of pregnancy) may result in abortion, stillbirth, or premature delivery with amnionitis (during late pregnancy). Neonatal infection can be fatal and may be a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 of early onset (occurring within the first 24 hours to first week of life) or late onset (7 days to 6 weeks after birth) and is usually nosocomially acquired. Early onset disease (EOD) accounts for 62% of cases and may manifest as sepsis, respiratory distress and thrombocytopenia, while late onset disease (LOD) may present as sepsis or meningitis. Case fatalities in neonates are high [4] (Okike et al., 2013).
The 2017/2018 outbreak of listeriosis in South Africa, highlighted the need for adequate surveillance and prevention of this infection, bearing in mind that prior to this outbreak, the NICD regarded Listeria monocytogenes as a rare pathogen in meningitis infection in South Africa (NICD 2017) [5].
Mathematical modeling has been utilised in the fight against various epidemics. In particular, there is a vast literature on so-called compartmental models of systems of ordinary differential equations (ODEs), [6][7][8][9] for instance. Such models are mainly utilized under the assumption that there are no external forces working on the system. Understanding exactly how a sudden epidemic outbreak occurs will often require a different analysis. Such analyses of outbreaks have been studied for various diseases, [10], [11] and [12].
The aim of this paper is to propose a compartmental model for the population dynamics of listeriosis, and to derive a model for analysing a listeriosis outbreak. In particular the model explicitly accommodates neonatal infections. Similarly as is commonly done for cholera models, for instance [7][8][9], we include a compartment to represent the reservoir of bacteria. We also include a compartment to represent the incubation phase. For the 2017/18 outbreak [13] that happened in South Africa, we calculate the time pattern and intensity of the force of infection, and we make projections on the numbers of new infections and mortalities subsequent to the recall of the contaminated food that was the source of the infection. The model is calibrated using data from the NICD [13] of South Africa, together with existing data in the open literature not necessarily from South Africa. Our model differs from both the listeriosis population dynamics model in [14] and the listeriosis-anthrax coinfection model in [6] in several ways. Essentially the latter two models do not explicitly consider neo-natal infections and do not include a compartment for the incubation phase.
There are different ordinary differential equations (ODE) models for the population dynamics of infectious diseases under different conditions. Usually such models assume that the infection has entered the population or the system, and thereafter there is no external interference on the system. We shall refer to such models simply as population models. Modelling of how the population responds (numerically) while the infection is invading the system, i.e., during an outbreak of the disease, requires a different approach. In this paper we propose a population model for listeriosis disease. Properly calibrated, the model permits computation of future outlooks on a listeriosis epidemic. In terms of this proposed model we reflect on the listeriosis outbreak in South Africa 2017-18 [13]. In particular we discover the temporal pattern description of how the original infection took place and we make future projections of its epidemiology. For a previous listeriosis outbreak, the paper [12] presents a similar investigation but with completely different methods.

Methods
The population dynamics of listeriosis disease is studied by way of a newly introduced compartmental model of ODE. The model explicitly accommodates neonatal infections and we include a compartment to represent the reservoir of bacteria. Recent examples of models with vertical transmission can be found in the papers [15] and [16] for instance. Our analysis includes both theoretical mathematics as well as computer simulations. The model is calibrated to South African data.

The model
We describe our ODE compartmental model as follows. By E we denote the collection of individuals who are infected with listeriosis, but do not experience the symptoms. Those who are beyond the latter phase and are suffering the symptoms of listeriosis form a compartment I. By R we denote the removed, i.e., those who have recovered from listeriosis and have temporary immunity. This can be bodily immunity or through lifestyle adaptation. The remaining population forms the class S of so-called susceptibles. At any given time t the sizes of these compartments are denoted by E(t), I(t), R(t) and S(t) respectively. The size of the population at time t is denoted by N(t), and then N(t) = S(t) + E(t) + I(t) + R(t). It is also assumed that the pathogen is present in the environment, over and above those that live in the form of human infection. Thus our system has another compartment B, with B(t) representing the abundance of the pathogen at time t. The data in the report [13] of the NICD show that in any given week the number of new infections were always below 1000, in a very large population. This has implications (mainly simplification) of our computations. Since we do not have data available to feed into the model, we work with B(t) just as an abstract quantity, but we now declare the (proportional) connection of B(t) with the concentration of the pathogen. If we write P(t) for the average (over the population) concentration of the pathogen in the diet, then for a suitable proportionality constant q, we write B(t) = qP(t). Here q is chosen in such a way that: B ¼ 1 corresponds to the value of P that; in a population of 60 million susceptibles; causes100 infections per week: ð1Þ The sizes of these compartments will be treated as continuous variables (rather than integers). We assume a homogeneously mixing population. We introduce non-negative constant parameters K, μ, ϕ, ψ, δ, z and ν. The number K is the size of the population when in diseasefree state, i.e., at disease-free equilibrium. The general mortality rate is denoted by μ. In the Icompartment the mortality rate is amplified by the disease-induced mortality rate δ. Therefore in the compartments S,E,R and I respectively, the rate of outflow due to death are (resp.) μS, μE, μR and (μ + δ)I.
It is assumed that there is an inflow of individuals (newborns) into the system at a rate μK. Part of the inflow is into the I-class, as neonates are at risk of contracting the infection from the mother if she is infected. So let us use θ to denote the fraction of neonates from infected mothers who do contract the infection from the mother. The number of babies from infected mothers can be calculated as μI. Therefore the number of newborns entering the I-class is θμI.
We assume that the probability of an individual becoming infected with listeriosis during a short period of for some fixed positive constants c and L. This means that new infections occur at a rate β(t)S (t), and so β(t)S(t) is the rate of transfer from the S-class to the E-class. The given form of β(t) is also applied in modeling of cholera. The rate of transfer from the E-class to the I-class is ϕE (t) and the rate of transfer from the I-class to the R-class is ψI(t). From the R-class, individuals move to the S-class at the rate zR(t). The flow chart, Fig 1, of the proposed model shows the rates of transfer in or out of compartments.
In the environment, the bacteria are assumed to have an effective rate of decline denoted by ν. Furthermore, the infectious class I of the human population contributes to the bacteria in the environment at a rate ξI(t). The factor ξ takes the form for some non-negative constants η and η 0 . Then the population dynamics of the disease is described by the system (3) of five ordinary differential equations.
dSðtÞ dt ¼ mK À ymIðtÞ À mSðtÞÞ À bðtÞSðtÞ þ zRðtÞ For notational convenience we introduce the symbol The state of the population at any time t is a solution X(t) = (S(t), E(t), I(t), R(t), B(t)) of this system of differential equations. As is common in compartmental ODE models, this system of equations is applicable when studying the phase subsequent to an outbreak. During the outbreak itself, the infection is imposed on susceptibles by an external force, and we follow a slightly different approach.
It can routinely be shown that the set is an invariant domain for the system. This means essentially that all the state variables will be non-negative and the total population size will never be bigger than K. The basic reproduction number The model has a disease-free equilibrium point, denoted by X � 0 , and X � 0 ¼ ðK; 0; 0; 0; 0Þ. The basic reproduction number R 0 of a disease in a population is defined as the number of infections caused by introducing a single infectious person into a susceptible population. It has the property that whenever R 0 < 1, then the disease-free equilibrium is locally asymptotically stable, see [17]. The latter means that if the initial state of the population is sufficiently close to disease-free, then in the long run the disease will be eradicated from the system. We determine the basic reproduction number of the model (3) following the next-generation matrix method as presented in the work [17]. The vector of infectious classes is (E, I, B). Then the next-generation matrix is of the form FV −1 with F and V as below. Therefore we obtain the next-generation matrix as: The basic reproduction number is the maximum of the moduli of the roots λ of the latter cubic equation. Therefore, ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ym 2m 1 Recall that the number R 0 serves as an indicator of local stability (or otherwise) of the diseasefree equilibrium of system (3) (see [17]). We shall also refer to the following number R g , when studying global stability, proving that the condition R g < 1 guarantees global stability of the disease-free equilibrium. Global stability of the disease-free equilibrium means that irrespective of the initial state of the population, in the long run the disease will tend towards extinction, and is a stronger condition than local stability.
In fact, as indicators of stability, the invariants R g and R 0 are equivalent as we prove below.
Proof. We can write ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ð1 À xÞ This proves (a).
(b) This follows by a similar argument.

Global stability of the disease-free equilibrium
When working towards eradication of an endemic infectious disease, global stabilization of the disease-free equilibrium state is important. This will ensure that in the long term, the pathogen will become extinct in the population. Global stability theorem. If R 0 < 1, then the disease free equilibrium, X � = (K, 0, 0, 0, 0), of the system (1) is globally asymptotically stable.
Proof. Let us assume that R 0 < 1. Then also R g < 1. This condition is equivalent to the inequality Given the latter inequality there exist positive numbers b 0 , b 1 , � such that We fix such numbers b 0 , b 1 , and �. Now we define the following numbers: In particular we note that ym À ðc þ m þ dÞ þ a 3 Z < 0: Therefore we can find a number a 0 > 0 such that a 0 � b 1 and such that ða 0 þ 1Þym À ðc þ m þ dÞ þ a 3 Z < 0: We can also choose a number a 2 > 0 such that We introduce a new variable Q(t) = K − S(t) and we define a function V ¼ VðQðtÞ; EðtÞ; IðtÞ; RðtÞ; BðtÞÞ; which we shall prove to be a Lyapunov function with respect to the disease-free equilibrium point X � 0 . This will suffice to complete the proof. Thus, let þ R½À a 2 ðz þ mÞ À a 0 z� À a 3 nB: Noting that β = cB/(L + B) and S � K, in particular we obtain the following inequality Also, ξ � η. Therefore we obtain: On the right hand side of the latter inequality, we note that the coefficients of Q, E and R are negative. From the inequality (5) it follows that also the coefficient of I is negative. Regarding the coefficient of B we observe: and so by the inequality (4) the coefficient of B is negative. Thus we have shown that V is a Lyapunov function, and this completes the proof.

Parameter values relevant to outbreak
Values of parameters of system (3) will be discussed in this section, and are presented in Table 1. We also analyse the transmission probabilities α(t) for the duration of the outbreak.
Since the listeriosis data [13] from the South African NICD are given per week, we choose 1 week as our unit of time (but in simulations we consistently use timestep equal to 1 day). From the NICD data which are reflected in Fig 3, we fit the parameters ϕ, ψ and δ, and their values are as in Table 1. The data also sheds light on θ. Regarding θ, simple calculations reveal that there must have been many infected babies whose mothers were not listeriosis patients. In a random sample of n the expected percentage F of babies, who in a given period of 62 weeks were less than 28 days old, is calculated as The percentage of neonates among the infected is reported in [13] as 42% which is bigger by a factor of more than 20. This leads one to conclude that if both a mother and her neonate was exposed to L. monocytogenes, then the baby was more likely to contract the disease than the mother. In particular then, for the given scenario one would suspect that if the mother of a neonate was infected, then her child would very likely also be infected. On this basis we guess a value, θ = 0.9. The other parameters are obtained from the literature (except for L, η and ν, which are not relevant to the outbreak).
In the census document [18] the average life expectancy in South Africa is given as 61.2 years for males and 66.7 years for females. Therefore we calculate the value of the general mortality rate μ for South Africa as 1 ð64:45Þð52Þ ¼ 0:000298 per week. Members in the class R are considered at least temporarily safe from becoming infected, due to bodily immunity or through lifestyle adaptation. In the case of listeriosis there are clear guidelines for avoiding infection, generally published regularly by public health authorities, see [13] or [20] for instance. The effective duration of the temporary immunity is taken as 6 weeks, and so z = 1/6 per week.
Regarding the function β, it is difficult to find an optimal value for L. However, we assign to it a nominal value L = 1. Then the value of c is fixed by the definition of B(t): Since B = 1 corresponds to 100 new latent infections per week and we have chosen L = 1, we obtain c × 1 × 6 × 10 7 /(1 + 1) = 100. This means that c = 1/300000 (per week).
The parameters η, η 0 and ν will be given nominal values. They are not relevant to the outbreak and we shall explain them at a later stage.
Remark (a) A value for ϕ can be estimated based on the papers [21] and [22]. The incubation period of listeriosis varies quite significantly: from under 24 hours up to as much as 90 days-a study of the incubation period appears in [21] for instance. A final value is not concluded in [21]. Our fitted value of ϕ, puts the average incubation period for listeriosis at just below 5.5 days.
(b) In the study [22] it is observed that the mortality among the infected is one-third. The period post incubation, subsequent to the pathogen having taken effect, over which these deaths occur was not specified. We note that treatment for listeriosis takes 2-6 weeks (Bacteremia should be treated for 2 weeks, meningitis for 3 weeks, and brain abscess for at least 6 weeks, [20]). So if one takes the duration of the infectious period as 6 weeks, then it suggests a disease-induced mortality rate of 1/(3 × 6) per week. This is approximately three times the rate that we calculate for this outbreak.

During the outbreak
In what follows we present a numerical analysis of the force of infection, α(t), during the outbreak period, 01 January 2017 up to 08 March 2018. We compare the model output of the mortalities and recoveries with the numbers reported in [13], and we present an outlook on the return of the system to equilibrium subsequently to the recall of the food products that caused the outbreak.
As noted earlier, the original system of equations is applicable when studying the phase subsequent to the outbreak, when there are no external forces on the system. During the outbreak the infection is imposed on susceptibles by an external force. This force is assumed to be resulting in infections at a rate α(t)S(t), where α(t) is a function of time. This α(t) replaces β(t) as the driving force of the infection, it does not have the same functional form as the β(t) of the population model and is unknown. In order to assess the force of infection, α(t), during the outbreak period itself, we shall utilise the system (7) below. System (7) is obtained from system (3) by excluding the fifth differential equations and the force of infection will be determined from the data [13]. Then we have a system (7) of four ordinary differential equations.
Sðt þ 1Þ ¼ SðtÞ þ ðmK À ymIðtÞ À mSðtÞÞ À bðtÞSðtÞ þ zRðtÞÞh Initial values for the compartmental variables are chosen as below, which we shall see eventually are the endemic equilibrium values corresponding to the parametrization given in Table 1: The new cases for each week can be read from the data in [13]. During the t'th week this number is the same as (ϕE(t) + θμI(t)), and the expression for E(t) contains the variable α. These can be used to calculate the function α(t), which we display in Fig 2. In order to smooth out the NICD data, we found a density plot, p(t), for the relative number of cases, i.e., (cases)/967. This density plot is shown together with the scatter diagram of relative case numbers in Fig 3. This p(t) will be used hereafter to describe the force of infection over the outbreak period, in a similar way as is done for the rough data.

Comparing the model output with reported results
In Table 2 we present three outputs of the model for comparison with reported results [13], for the period 01 January 2017 up to 08 March 2018. We give the total new cases as reproduced on the basis of the force of infection, computed using the Eq (8) above. The three outputs that we compute are the total number of cases, the total number of fatalities and the total number of recoveries. We are then able to calculate the number of remaining infectives through simple arithmetic. The outputs are calculated twice: i.e., first for rough data and then from the density plot. In both cases the outcomes compare very well with the observed data.
In Fig 4 (based on the NICD-data) and Fig 5 (based on the density plot) we show the trajectories of the E-, I-, and R-compartments over the outbreak period. The S-class is many orders bigger. In fact, we note that the total number of cases during the outbreak was fewer than 1000. Therefore if we consider S(0) to be 60 million, then during the outbreak S(t) will always be above 59999000.   [19] announced that NICD records show between 60 and 80 mortalities p.a. due to listeriosis, over the years 2012-2016. Let us take this as approximately 70 listeriosis deaths p.a. under equilibrium conditions. This enables us to calculate I � as follows: δI � × 52 = 70, therefore I � = 71.6. We note that at endemic equilibrium, from the second and the third equations of the system (3) we have:

PLOS ONE
This yields Therefore we have calculated the value of B � in terms of the known parameters, the values of which are well motivated. In particular we did not require ν, η or η 0 . Thus we are in a position to calculate the ratio ξ � /ν = B � /I � . This yields ξ � /ν = 0.000342. So, in order to complete the calibration of the model, we only need to determine the value of either ξ � or ν. For both of the latter two parameters, their values depend on the environment together with human behaviour and lifestyle. So for any two different populations, there may be a huge difference between the two η -values. Likewise for ν -values. We pick a nominal value ν = 2 units per week. Then we use η 0 = 1/71.6 and that automatically determines η: Now we can harness the model for making future projections. The graphs in  Table 3 shows the equilibrium values as calculated using the formulae and parameter values, and also the values of the variables S, E and I as predicted by the model for the dates 30 June

Reducing the annual mortalities
The stability threshold is computed as having the numerical value R 0 = 1.43. The Global Stability Theorem revealed that when R 0 < 1, then in the long run, the disease will vanish from the system. This informs how much should be done by way of intervention to ensure that the disease will not persist and flourish. Interventions should lead to reduction of R 0 below unity.
In particular, with an intervention which changes the parameter values in such a way as to result in R 0 going below unity, the disease will ultimately vanish from the population. In this regard we could immediately consider parameters that depend on behaviour or lifestyle, such as η: perhaps through improved sanitation service, or c: good hygiene around food and drinking water. Indeed the graph in Fig 7 illustrates, how the infectious class diminishes converging to zero, under a suitable intervention. In this case, the parameter values are all the same as in   Table 1, except for η which has been reduced to 0.00065 (instead of 0.0014), and this gives R 0 = 0.987 which is below 1.

Discussion and conclusion
We have constructed a mathematical model for the population dynamics of listeriosis infection. Such a model is useful to public health authorities when exploring intervention strategies to control or eradicate an infectious disease (Huppert and Katriel 2013 [23], Althaus et al 2014 [24]; Rose et al 2014 [25], Heesterbeek et al 2015 [26]). In particular, interventions should aim at reducing the value of the basic reproduction number, bringing it to below unity where possible. In the case of the current model, having R 0 < 1 will guarantee eradication of the disease in the long run. For some disease models, the said condition is not sufficient to ensure eradication, but we proved that for the current model it is indeed. We were able to harness the model in a modified form, retrospectively to reveal the intensity, or at least the temporal pattern, of the infection which caused the outbreak. The data available is limited, and consequently not all of the parameters could be determined. Nevertheless we were able to obtain a reasonable estimate of the basic reproduction number and the endemic equilibrium point. The endemic equilibrium as calculated yields an endemic equilibrium value of I � which is in line with the annual infection numbers in South Africa (see [19]). We also demonstrated how to compute future projections (Table 3), and how to test an intervention for the purpose of eradication (Fig 7). With a change in the value of η to η = 0.00065, we obtain R 0 = 0.987 which is below 1, and as guaranteed by the Global Stability Theorem, over time the infection will be eradicated. The graph shows how the system will change over a period of 130 weeks, starting from what is currently the endemic equilibrium. We note a definitive tendency of convergence to disease-free state in the longer term. For application to another population, the model can be used in the same form, but the parameter values must be adapted. For this reason adequate surveillance systems and resources to diagnose or detect listeria must be in place. Accuracy can be improved by allowing for age structure and spatiality, and to produce simulations through more accurate methods. Climate variability may directly or indirectly lead to higher infection rates of listeriosis in South Africa [27]. Modeling that accommodates the effect of climatic factors, such as was done for malaria by several authors ( [28] or [29] for instance) may prove useful in the near future.