Fundamentally, infectious diseases persist in animal populations because the normal patterns of social contact required for day-to-day life create repeated opportunities for transmission from infected to susceptible individuals. By calculating a value for R0, we can estimate the likelihood that an infectious agent will establish and persist in a host population. In the previous module, we represented transmission as a continuous-time process, where individuals transition between disease states at small rates each day. For populations such as continuous-flow swine breeding systems, that can be a reasonable approximation. However, in other systems such as pastoral sheep and beef breeding herds, key demographic processes (births, aging, pregnancy, and removals) tend to occur in seasonal pulses rather than continuously. In extensive pasture-based systems, animals may also be less well-mixed, with contacts shaped by paddock structure, social grouping, and management practices, which alters within-herd transmission dynamics. In this module, we will explore how to adapt our models to reflect seasonal demography and more realistic patterns of contact within herds.
By the end of this module, students should be able to:
In the previous models developed for Module 2, we assumed simply that animals within the population were all homogenous and had equal probabilities of acquiring and transmitting disease. The situation gets much more complicated for situations when disease transmission occurs at different rates or where the clinical outcomes of the disease vary depending on the demographic characteristics, management groupings, and health status of different individuals within a population. For example, in the case of seasonal influenza for humans, we know that children are much more infectious than adults and so we would want to develop a model that tracks these two demographic groups separately. Fortunately, there are relatively straightforward techniques we can use to modify our compartmental models to account for these dynamics.
In the above system, note that we have (1) created a separate set of SIR compartments for children and adults, (2) modelled the infection pressure for each susceptible compartment as a function of both the pressure from children (βchildIchild) and the pressure from adults (βadultIadult) with different transmission coefficients for each, and (3) described the aging process with vertical arrows linking each SIR compartment for children with the corresponding SIR compartment for adults and described the rate of aging as ρ. Just make sure that the aging rate corresponds with the time step of your model (i.e. if your model is updated on a time step of days and you classify a child as being under 10 years then the rate ρ would be (1/3650 days) and not (1/10 years).)
Read through the description of bovine viral diarrhoea virus pathogenesis below:
Bovine viral diarrhoea virus (BVDV), the causal agent of BVD and mucosal disease complex, is classified in the genus Pestivirus in the family Flaviviridae. Disease induced by BVDV varies in severity, duration, and organ systems involved. Infection of immunocompetent susceptible animals with either noncytopathic or cytopathic BVDV is termed acute or transient BVD. Inapparent or subclinical infection without any clinical signs that is followed by seroconversion is the most common form of infection in the field. Acute clinical disease may range from mild disease of high morbidity and low mortality to severe enteric disease with considerable mortality. Biphasic fever (~104°F [40°C]), depression, decreased milk production, transient inappetence, rapid respiration, excessive nasal secretion, excessive lacrimation, and diarrhoea are typical signs of acute clinical BVD. Clinical signs of disease usually are seen 6–12 days after infection and last 1–3 days. Transient leukopenia may be seen with onset of signs of disease. Recovery is rapid and coincides with production of viral neutralizing antibody that confers lifelong immunity against the disease. Gross lesions seldom are seen in cases of mild disease. Lymphoid tissue is a primary target for replication of BVDV, which may lead to immunosuppression and enhanced severity of intercurrent infections.
In pregnant cattle, BVDV may cross the placental barrier and infect the fetus. The consequences of foetal infection usually are seen several weeks to months after infection of the dam and depend on the stage of foetal development and on the strain of BVDV. Persistent infection develops when noncytopathic BVDV is transmitted transplacentally during the first 4 mo of foetal development. The calf is born infected with virus, remains infected for life, and usually is immunotolerant to the resident noncytopathic virus. Transplacental infection that occurs later in gestation can result in abortion, congenital malformations, or birth of healthy calves that have antibody against BVDV. The prevalence of persistent infection varies among countries and between regions within a country. The estimated mean animal prevalence of persistent infection with BVDV is ~1%–2% but may approach 4% on dairy farms with endemic BVDV infection. On a given farm, persistently infected cattle are often found in cohorts of animals that are approximately the same age. Persistently infected cattle can shed large amounts of BVDV in their secretions and excretions and readily transmit virus to susceptible herdmate compared with transiently infected animals that shed only low amounts of the virus.
Calves born to recovered dams will general have antibodies against BVD for 6 to 10 months. Calves born to persistently infected dams will always be persistently infected themselves.
| Symbol | Definition |
|---|---|
| κ | Rate of loss of maternal immunity |
| βPI | Transmission from PI animals |
| βTI | Transmission from TI animals |
| γ | Recovery rate |
| α | Increased mortality rate in PIs |
The effects of BVD on animal health and production depend on the animals stage of gestation at the time of infection so we might want to include demographic compartments describing early pregnant cows, mid pregnant cows, late pregnant cows, and open cows.
Depending on the herd management structure, you might also want to include separate compartments for calves (0 to 6 months of age), replacement heifers (6 to 15 months of age), and breeding cows (15 months to cull) since the transmission rate within a management group is likely to be much higher than transmission between groups.
You would need to create separate disease compartments for each demographic group and describe the transition between the demographic states (vertical arrows) as well as the transition between the disease states (horizontal arrows). When it starts to get this complex, I will often switch to an individual based model where you just have to track disease status and demographic status for each animal.
For many decades, R0 has been considered the holy grail of measurements in infectious disease epidemiology because it gives an indication of spreading potential of pathogen and can be used, for example, to estimate the percentage of the population that needs to be vaccinated in order to achieve herd immunity. There are five main ways we can generate estimates for R0:
If we have outbreak data showing the number of infected individuals over time, then we can approximate R0 from the initial multiplication rate using the following formula:
Where D is the length of the infectious period and td is the initial doubling time. This assumes that the population of susceptible animals is unlimited.
If the proportion of susceptible individuals remaining in the population after an outbreak is known, we can estimate Ro as:
where Se is the total number of susceptible individuals at endemic equilibrium and N is the total population size. For example, if 100% were initially susceptible and 20% remained susceptible until after the outbreak, the approximation yields Ro = 1/0.2 = 5. If we were trying to estimate this from real data, we could measure the seroprevalence of disease in the population and assume that everyone else is susceptible (i.e. if 60% of the population have antibodies against the disease, Ro = 1 / (1 – 0.6) = 1.67). This method of calculation is applicable when populations are dynamic and new susceptible animals are added as mortality reduces the population. If an equilibrium is not established, the proportion of susceptible animals left after the outbreak is under-estimated and Ro becomes over-estimated.
A formula for a situation where not all individuals were initially susceptible is:
Where μ0 is the proportion of individuals that were susceptible at the start of the outbreak and μα is the proportion of individuals that were susceptible at the end.
If infection is strongly related to age and typically persists for life, Ro may be estimated by demographic data. For example, calves are infected with M. paratuberculosis (Johne’s disease) between birth and about 12 months of age and remain infected for life (even though the age at shedding may start only years after infection). Under such circumstances, Ro can be calculated as follows:
Type I: Ro = L/A
Type II: Ro = 1 + L/A
Where A is the average age at infection and L is expected lifetime. The types I and II refer to the type of survival curve for the age class of animals: type I is a fixed lifetime with constant survival while type II is an exponential survival curve. Up to 2.5 years of age, the death rate of Johne’s-infected cattle is very low, hence the rectangular type I survival curve, while selective culling after this age leads to an exponential decrease of type II survival.
Deriving an equation for R0 is relatively straightforward when we have simple infection dynamics such as the endemic SIR model below.
We estimate R0 simply by dividing the rate individuals get infected by the rate individuals are recovered or removed. Remember that at time of initial disease invasion, all individuals in the population (N) will be susceptible.
This makes biological sense because in order for the disease to persist the rate that individuals are getting infected must be greater than the rate that they are recovering.
Things get a little bit more complicated when we are dealing with diseases that have multiple compartments containing infected individuals. For example, this commonly occurs when there are latent stages (E) or carrier stages (C) in the model infection dynamics or when we are dealing with populations that are structured by age, sex, species, or behavioural risk. The number infected individuals at any given time depends on the rate that susceptible individuals initially become infected (transmission) as well as the rate that infected individuals move through the different infection states (transition).
Fortunately (or unfortunately depending on how much you like matrix algebra), there are other methods we can use to derive estimates of R0 for complicated models. If you struggle to understand matrix algebra like I do, there are many user-friendly tutorials on the Khan Academy website:
Our example will be an SEIR model that includes births and deaths. As shown below, there are four compartments in the model and two of these compartments are infected.
In the SEIR model, individuals can enter the E compartment as new infections (transmission event), move from the E compartment to the I compartment (transition event), and transfer out of the infected compartments (E and I) via natural death and/or recovery (exit events). Only infectious (I) individuals are capable of transmitting disease even though both E individuals and I individuals carry the pathogen.
For any infected compartment, the rate of change in the number of individuals in the compartment is given by the formula:
where x denotes the number of individuals in compartment i and :
We can rewrite the expression as:
We call Fi(x) the transmission matrix that describes how naïve individuals become infected and Vi(x) the transition matrix that describes how infected individuals move between and out of infected states. For various mathematical reasons (which I freely admit to not understanding at all) when you multiply the transmission matrix (F) by the inverse of the transition matrix (V-1), the dominant eigenvalue of the resulting matrix becomes the Ro equation. Don’t panic – it sounds really complicated, but the steps to generate the Ro equation are actually straightforward to follow.
Remember that we are only interested in the infected compartments at this stage. Since there are two of them in our system, we will need to create a 2×2 matrix and fill in the appropriate rates (aij). You can think of the subscript notation aEE as the rate a susceptible individual moves into the E compartment as the result of contact with an E individual, aEI as the rate a susceptible individual moves into the E compartment as the result of contact with an I individual, and so on.
The only new infections in this case are entries into compartment E as the result of an infectious contact between an individual in compartment S and an individual in compartment I (aEI), which occurs at the rate βIS.
Remember that at the disease-free equilibrium the total number of susceptible S individuals is equal to N and the number of infectious individuals is equal to 1. Making that substitution, we get the complete matrix F:
We next follow almost exactly the same procedure to generate the transition matrix, which has two types of events: transition events and exit events.
In this case, you can think of the subscript notations aEE as an individual moving out of the E compartment, aEI as an individual moving into the E compartment from the I compartment, aII as an individual moving out of the I compartment, aIE as an individual moving into the I compartment from the E compartment.
There are no transitions into E (remember that βSI represents a transmission event) and individuals leave E through either death (μ) or transition to I (δ). Individuals transition into I from E as they become actively infectious (δ) and then leave I through death (μ) or recovery (α). Making that substitution, we get the complete matrix V:
Note that the top-left and bottom-right cells are always positive, while the top-right and bottom-left cells are always negative. We also don’t need to worry about the multiplying by the compartment letters (this is because there is another underlying step of taking the partial derivatives of the terms in the matrix).
Now that we have defined the two matrices F and V, we are ready to create the next-generation matrix FV-1. This involves a lot of complicated maths steps. If you are interested in reading more technical and mathematical explanations of the next-generation matrix method for calculating R0, check out the following references:
Fortunately, there is a R package (‘rSymPy’) that has some easy functions to help us with the underlying math, which are outlined in the script file (R – Next Generation Method Ro). Note that the package will only run in the 32-bit version of R due to issues with how it handles another required Java library.
This produces the following output:
and all you have to do is re-substitute the symbols for the letters to generate your Ro equation.
As you can see from the equation, R0 becomes a function of the rate that individuals enter and exit the infected E and I compartments. If we have parameter estimates for β, µ, α, and δ, we can plug them into the equation to get the numerical point estimate for R0.
For the next example, we will look at an SEIRC model that has three infected compartments: exposed (E), infectious (I), and carrier (C). We will assume that the transmission rate β is the same for infectious and carrier individuals. Since there are three infected compartments, we will need to construct 3×3 matrices for F and V.
There were only two transmission events in this example: S individuals becoming E as the result of contact with I individuals and S individuals becoming E as a result of contact with C individuals.
We will next create our new 3×3 matrix for V and then populate the values based on our model equations.
Now that we have our F and V matrices defined, we will need to load the rSymPy package in R.
> install.packages("rSymPy")
> library("rSymPy")
This package provides tools that allow you to substitute symbols for numeric values in common mathematical functions. This is referred to as symbolic algebra or symbolic computation. NB – the syntax will look different than you are used to in R because the package interfaces with another programming language (Python) to run the calculations.
To make the coding easier, we will temporarily assign the elements in our F and V matrices with unique single letter symbols.
We first need to let rSymPy know what symbols are used in the matrices.
> sympy("var('a,b,c,d,e,f,g')")
Then we define our recoded F and V matrices. Note that rSymPy syntax requires you to enter the matrix data row by row. The outer cat() function is just to make the R-console output pretty.
> cat(sympy("F = Matrix([[0,a,b], [0,0,0], [0,0,0]])"), "\n")
> cat(sympy("V = Matrix([[c,0,0], [-d,e,0], [0,-f,g]])"), "\n")
We can now multiple F by V-1 to get the next-generation matrix.
> cat(sympy("N = F*V.inv()"), "\n")
[a*d/(c*e) + b*d*f/(c*e*g), a/e + b*f/(e*g), b/g]
[ 0, 0, 0]
[ 0, 0, 0]
Conveniently, there is also a function that will produce a vector of eigenvalues for the matrix. This will allow us to easily identify and select the largest eigenvalue as our equation for R0. The output below means that the eigenvalue a*d/(c*e) + b*d*f/(c*e*g) occurs once while the eigenvalue 0 occurs twice. Remember there are 3 eigenvalues for a 3×3 matrix. Obviously, in this case, the first one is our dominant eigenvalue.
> sympy("N.eigenvals()")
[1] "{a*d/(c*e) + b*d*f/(c*e*g): 1, 0: 2}"
The last step is to go back and substitute the real matrix elements for the letter symbols. This makes the final equation for R0:
Glad we didn’t try to work this one out by hand! We can again substitute the values for β, µ, α, θ, ρ, and δ to get the point estimate for R0.
In compartmental models, demographic processes have traditionally been modelled as continuous processes or rates where fractions of animals move between compartments every time step. However, we know that management events such as culling often occur at a single time point (discrete processes) or, in seasonal herds, animals will only be pregnant during certain times of year (seasonality). It is possible to write complex equations for the transmission coefficient that will allow the demographic and disease dynamics to fluctuate over time while still keeping the model as a series of equations that can be solved.
However, we can also actually very easily build discrete management events into a traditional compartmental model – the only problem with doing this is that we can no longer use some of the standard methods for estimating a value of R0 from the model equations since we are no longer modelling the entire system as a single set of equations. That’s not necessarily a problem because the concept of R0 is not straightforward to apply to industry level metapopulation models where is there is no such thing as an “average” individual interacting with a completely susceptible population.
Advancing Animal Care through Shared Learning