|Home | About | Journals | Submit | Contact Us | Français|
This paper investigates the early viral dynamics of foot-and-mouth disease (FMD) within infected pigs. Using an existing within-host model, we investigate whether individual variation can be explained by the effect of the initial dose of FMD virus. To do this, we consider the experimental data on the concentration of FMD virus genomes in the blood (viral load). In this experiment, 12 pigs were inoculated with one of three different doses of FMD virus: low; medium; or high. Measurements of the viral load were recorded over a time course of approximately 11 days for every 8 hours. The model is a set of deterministic differential equations with the following variables: viral load; virus in the interstitial space; and the proportion of epithelial cells available for infection, infected and uninfected. The model was fitted to the data for each animal individually and also simultaneously over all animals varying only the initial dose. We show that the general trend in the data can be explained by varying only the initial dose. The higher the initial dose the earlier the development of a detectable viral load.
Foot-and-mouth disease (FMD) is an infectious viral disease (genus: Aphthovirus, family: Picornaviridae) that affects cloven-hoofed animals, such as cattle, pigs, sheep and goats. It has an almost worldwide distribution with an endemic and occasionally epidemic occurrence. If an outbreak occurs, it can have a significant economic impact on the areas affected. For example, in the UK, the 2001 FMD outbreak was estimated to have caused a loss of £3.1 billion to the agricultural sector and food chain, with a similar amount lost by the tourism industry (Thompson et al. 2002).
The spread of FMD has already been modelled at different levels, in particular farm-to-farm transmissions and within-farm transmissions. A general discussion of these types of models can be found in Woolhouse (2004). One area of FMD dynamics that has received little attention is a mathematical description of the viral dynamics within a host. An understanding of the within-host viral dynamics can be used to aid the study of FMD spreading within a population, providing more accurate information, such as incubation and infectious periods. Quantitative descriptions of FMD viraemia have already been reported (Cottral & Bachrach 1968; Sutmoller & McVicar 1976; Alexandersen et al. 2001, 2003b; Quan et al. 2004), as well as the clearance of virus from the circulation (Sutmoller & McVicar 1976).
Within-host modelling of viral dynamics has previously been directed at human diseases, in particular human immunodeficiency virus (Perelson et al. 1996, 1997; Corfec et al. 2000; Nowak & May 2000), hepatitis B virus (Herz et al. 1996; Nowak et al. 1997; Nowak & May 2000; Lewin et al. 2001; Perelson 2002; Ciupe et al. 2007), hepatitis C virus (Perelson et al. 2005), influenza A virus (Baccam et al. 2006) and human T-cell leukaemia virus (Wodarz et al. 1999). These models have mostly focused on the chronic nature of the diseases, as well as chemotherapeutic interventions. In the models, the life cycle of the viruses is restricted to the circulation or closely associated organs, and there are no divisions between the circulation and organs. There have been no reports on the interaction of virus in the circulation and virus replication within organs. By contrast, pharmacokinetic models of the movement of drugs between different compartments within the body are well established (Riviere 1999).
In this paper, we examine the relationship between the initial dose of FMD virus (FMDV) injected into the blood and the subsequent within-host dynamics. In particular, we question whether individual variation is mostly due to the effect of the initial dose of FMDV. The basis of this study is data from an experiment where 12 pigs were inoculated with one of three different doses of FMD virus: low; medium; or high. Readings of the viral load were recorded over a time course of up to 11 days for every 8 hours. Our model then adds to the simple inspection of, and the statistical analysis of, the raw data in three important ways. First, the model suggests a possible mechanistic explanation for the observed behaviour. Second, it suggests that the wide range of dynamics seen in the experimental data can be largely explained by the differences solely in the initial dose. Third, the model makes testable predictions that can be investigated in subsequent experiments.
The model which we have chosen to use for this study was developed earlier (Quan 2005) to reflect the current knowledge of FMDV biology. We do not consider here alternative and more speculative models of the within-host dynamics, but show that this simple model can be used to explain the patterns shown in the experimental data by varying only the initial dose while also providing estimates of within-host parameters.
Twelve female white Landrace pigs, weighing 20–30kg, were split into three groups of four animals. The first group was housed in two separate boxes, each box containing two animals that were physically separated from each other by double, solid, watertight, wooden partitions (30cm apart), so that direct contact between the animals was not possible (Alexandersen et al. 2002a,b). The group was inoculated with a low dose of FMDV (106.1 FMDV genomes).
The second group was treated in the same way as the first except for receiving a medium inoculation (107.1 FMDV genomes).
The third group received a high inoculation (108.1 FMDV genomes), and the pigs were all housed in the same box without any separation between animals. This was due to the limitation of space in the isolation unit. Each animal had serum collected every 8 hours for up to 11 days.
The animals were housed in a secure isolation unit located at the Institute for Animal Health, Pirbright Laboratory. The unit consisted of five boxes in a row connected by a ‘clean’ corridor on one side and a ‘dirty’ corridor on the other. The personnel entrance to each box was from the clean corridor through the double doors with a shower unit between the doors. Before leaving the box, the personnel underwent disinfection with FAM (Evans Vanodine International, Preston, UK), followed by showering. Each box was self-contained with its own high-efficiency particulate air filtration, so that cross contamination was not possible.
The animals were inoculated intravenously in the anterior vena cava/brachiocephalic vein with 0.5ml diluted inoculum, either 1 in 10, 1 in 100 or 1 in 1000. This was equal to 108.1, 107.1 or 106.1 FMDV genomes, respectively. The inoculum used was isolated from a pig with FMDV strain O UKG 34/01 at Cheale Meats Abattoir, Brentwood, Essex, on 20 February 2001 (Alexandersen & Donaldson 2002; Alexandersen et al. 2002a,b, 2003a,c). The inoculum contained 107.2 TCID50 in the primary bovine thyroid (BTY) cells or 109.4 FMDV genomesml−1, and was diluted in HEPES tissue culture medium just before use. After inoculation, the dilutions were tested for genome quantity and infectivity using fluorogenic reverse transcription polymerase chain reaction (RT-PCR) and viral titration in the BTY cells.
Serum was collected by venipuncture (anterior vena cava/brachiocephalic vein) into plain serum tubes with no additive (Vacutainer; Becton Dickinson, Plymouth, UK). Samples were taken before (to act as controls) and then approximately 5min after inoculation to estimate the true nature of the initial dose of FMDV. Thereafter, the samples were taken every 8 hours: between 7.00 and 08.00, 15.00 and 16.00, and 23.00 and 24.00.
The virus was taken to be well mixed within the blood after 5min. Taking the blood volume to be 1/13 of the pig weight (Chang et al. 2002) gives a volume of 1.923l for a 25kg pig. The mean cardiac output of a typical pig is 4.5lmin−1 (Åberg et al. 2004). This gives an estimate of (1.923/4.5)×60=25.6s for the blood to circulate.
A detailed description of RNA isolation, reverse transcription, PCR and quantitation of the serum samples from the experiment can be found in the previous work (Quan et al. 2004).
The main assumptions in the model are as follows.
The model is a set of differential equations that describe the interactions between the virus concentration in the blood, X, and interstitial space, Y, and the proportion of available epithelial cells for infection that are either infected, F, or uninfected, U. The proportion of cells varies from 0 to 1; to obtain the number of cells either infected or not, the variables U and F must be multiplied by the total number of epithelial cells available for infection. A schematic of the model is shown in figure 1 and the variables X, Y, U and F are defined in table 1. The model differential equations are
where the initial condition X(0) is a fitted parameter, Y(0)=0, U(0)=1 and F(0)=0.
The ratio of volume of the blood to interstitial space, parameter r, was calculated from the volume of plasma (6.91% body weight) to interstitial space (13.09% body weight) (Pond & Hought 1978). Therefore, parameter r was fixed at a value of 0.53 in all our fits.
The lifespan of an infected cell, lf−1, was calculated as the approximate time for cell attachment and intracellular infection time of FMDV between 2.5 and 3.5 hours (Cartwright et al. 1957). For our fits, we chose to fix lf at 0.33, so that lf−1 is approximately 3.03 hours.
Initially, the proportion of infected cells, F(0), is zero. The initial concentration of the FMDV genomes in the interstitial space, Y(0), is also zero. The initial concentration of the FMDV genomes in the blood, X(0), is a non-zero value given by intravenous inoculation, and is treated as a parameter due to the variable nature of the procedure. The virus then spreads as follows.
The modelled dynamics are illustrated in figure 2. Sensitivity analysis of each parameter was investigated by scaling the parameter with different factors and plotting the resultant graphs.
The system of differential equations was solved using a fourth-order Runge–Kutta method. An adaptive step-size control was used to ensure the accuracy of the solution, which was necessary due to the rapid changes in the derivative of the solution (Press et al. 1986).
Parameters were estimated using maximum-likelihood methods, assuming lognormal errors in the data and left-censoring observations below a minimum detection threshold of 102 genomesml−1. The likelihood for the data, D, and parameter set, P, is given by
where Di is the observed concentration; xi is the expected concentration given by the model; and and Φ are the probability and cumulative density functions for a normal distribution, respectively. The number of data points is given by n. The observations are left censored as indicated by ci, where ci is 0 if Di>102 and 1 otherwise. The standard deviation of the measurement error, σ, was unfortunately not available, and so was estimated by considering it as a parameter that must be fitted (Stewart et al. 1998).
The likelihood surface was sampled using Markov chain Monte Carlo with 200000 iterations after a 200000-iteration burn-in period (see Gilks et al. 1995). Proposal distributions were given by normal distributions, where the standard deviations were adjusted before the burn-in period to achieve good acceptance rates using a Metropolis–Hastings acceptance rule. The sample with the highest likelihood was taken as the estimate of the maximum likelihood. The Nelder–Mead downhill simplex method (Press et al. 1986) was used to search the parameter space to find suitable regions in the parameter space to initialize the chain. The 95% CI of a parameter was taken as the range between the 2.5 and 97.5 percentiles of the estimated cumulative distribution function of that parameter. The model parameter analysis was calculated on a Unix platform using GNU-complied C++ code.
The data from each animal were fitted to the model individually. The model was also fitted simultaneously to all the 12 animals, where the parameters are constant across all animals, except for the initial dose of FMDV inoculated, X(0), which was fitted to each animal individually.
The likelihood for the common fit, com, is therefore given by
where Di and Xi are the data and the initial dose, Xi(0), for animal i, respectively, and P′ denotes the model parameters without the initial dose.
In our study, we wished to explore the early within-host dynamics of FMDV in pigs using a within-host model. We can make some observations of the data, as shown in table 2, before considering the model.
The reading recorded at 0.1 hours gives an indication of the effective dose that each animal received, since this was taken shortly after inoculation. First, consider the low-dose animals, UQ10–UQ13. Animal UQ11 had a smaller reading at time 0.1 hours, which subsequently has a late rise in the viral load. Second, in the medium-dose animals, UQ14–UQ17, animal UQ17 has a lower reading at time 0.1 hours than any other animals, which again shows a late rise in the viral load. Interestingly, this initial value is still higher than all the readings from the low-dose animals. The high-dose animals, UQ18–UQ21, all have high readings at time 0.1 hours.
The early, smaller peak in the data for the animal UQ15 indicates an early round of viral replication and is frequently seen in naturally infected animals. In this case, FMDV replicates locally before disseminating systematically. In this animal, it is possible that epithelial cells were infected directly when the needle was pushed and pulled out of the skin. However, we do not attempt to model this feature with our simple model, which aims to still capture the basic dynamics.
The viral load readings could not be reliably detected for the values below a threshold of 102 FMDV genomesml−1. This creates some uncertainty in the data, in particular with animals UQ11 and UQ17. The viral load readings for these animals are mostly low; therefore, the dynamics of the virus reproduction are mostly unobservable. The best-fit parameters in terms of the maximum likelihood for UQ11 are perhaps unrealistic, having a large delay in effective virus production and an ineffective clearance of virus after the peak. This fit is due to two factors in the data: first, the large number of values below the threshold of 102, including the initial value; second, the high later values with no clear peak. The model fit of UQ17 shows a delayed peak, which is likely to be the true behaviour of UQ11, which is supported by the model fitted simultaneously to all datasets.
Unfortunately, there are no data over the full time course of the experiment for every animal. Many animals were killed or died before the end of the experiment. Some animals clearly reached peak viral load: UQ10; UQ14; and UQ15. The data for other animals suggest that the peak viral load may have been reached: UQ13; UQ16; UQ20; and UQ21. The results for the remaining animals are more uncertain.
The best-fit parameter values for each animal are shown in table 3 and the corresponding 95% CIs in figures 5 and and6,6, calculated as described in §2.9. The resultant viral load curves of the best-fit parameters for each animal are illustrated in figure 3.
The model was fitted to all the 12 animals simultaneously, where the initial dose of FMDV was fitted to each animal, as described in §2.9. The common maximum-likelihood parameter estimates found are listed in table 3.
The fitted initial values of X for each animal are shown in figure 5b. The 95% CIs are illustrated with the individually fitted parameters in figure 6 (rightmost value). The viral load curves for the common fit are shown in figure 4. Note that the general shape of the curves is identical only to the timing of the peak differing.
Although there will be differences between the animals (such as weight), the model fitted over all animals simultaneously shows that the overall trend can be represented quite accurately by only varying the initial dose of FMDV (figure 4). Some compromise must be made between the individual model fits in order to obtain a good overall fit to all datasets. Therefore, the common fit of the model for each animal will not be as accurate as the model fitted individually. Nevertheless, the common fit captures the basic dynamics of each animal correctly, which is discussed further within §3.2.
The two model-fitting approaches were compared using Akaike's ‘An Information Criterion’ with the second-order information correction (AICc) (Sugiura 1978; Hurvich & Chih-Ling 1989, 1995), showing that the common fit model is the better approach, because it is using only 20 parameters compared with 96 parameters of the individual fits. This test accounts for the large difference in the number of parameters, as well as considering the small amount of data. (The individual fits have an AICc value of 825.7 and the common fit a value of 640.1.)
In this section, we consider each model parameter in turn: the role they have in the model dynamics and how they differ when fitted to different animals. The manner in which the model dynamics alter when a parameter is scaled is shown in figure 7.
If the flow of virus is decreased from the blood (kxy) to the interstitial space then few viruses are able to infect the cells creating a delay in the virus production. As more viruses are retained in the blood, an increase in the virus in the blood is also observed (figure 7a). Increasing the virus flow in the other direction (kyx) increases the virus in the blood since few viruses are cleared via the interstitial space. This then increases the infection rate of the cells as the virus re-enters the interstitial space (figure 7b).
The model may be simplified since the variables for the blood, X, and the interstitial space, Y, are in equilibrium with one another, related by the equation , although the results presented here refer to the full model (see appendix A). The sensitivity of changing the ratio between the blood and the interstitial space, given by constant r, is shown in figure 7c. Parameters kxy and kyx, which are fitted, can account for any variability in r. Therefore, we are not too concerned with the individual values of kxy, kyx and r.
The parameter ly controls the exit of virus from the interstitial space and therefore also controls the exit of virus from the blood. This parameter has two main roles: first to create the initial dip in the viral load, and second to create a decrease in the viral load after the peak.
As the value of ly is increased, the virus in the interstitial space is cleared quicker delaying the production of virus within the proportion of cells. The magnitude of the peak of the viral load also decreases since the rate at which virus is produced within the proportion of cells is met sooner by the clearance rate (figure 7d). For all of the individual fits, except UQ11 and UQ19, the upper bound of the 95% CI is unbounded. This is due to the exit rate to the blood, kyx, also being a variable and can be set as high as required, so that all of the viruses are not lost through clearance. Unfortunately, this means that only a lower bound is found for ly for these animals, showing that an almost instant clearance of virus from the interstitial space cannot be ruled out. The exit rate for the model fit of UQ11 is virtually zero, shown by the flat peak. The model fit of UQ19 relies less on virus production via the interstitial space and is therefore limited in magnitude.
The parameter hy controls the rate at which the proportion of cells becomes infected by the free virion population in the interstitial space. As the value of hy is decreased the production of virus is delayed (figure 7e), and if it were set to zero the virus would quickly disappear. Since the virus produced from the proportion of cells is returned to the interstitial space, this creates a positive feedback loop and the virus increases exponentially in the proportion of cells and the interstitial space. Lower values of hy create a larger initial dip and delay the onset of virus production (as mentioned above for the animal UQ14).
We can calculate typical values for the interstitial space infection rate of the proportion of available epithelial cells for infection, given by the expression hy(Y2/(m2 + Y2))YU, using the common fit model 24 hours before the peak in the viral load, at the peak and 24 hours after the peak, as shown in table 3 (animal UQ10 was used with a peak time of 138.64 hours). Since these values are for the proportion of available epithelial cells for infection, these values must be multiplied by the total number of available epithelial cells for infection to obtain values of the number of individual cells infected per hour.
The parameter hf controls the rate at which the proportion of infected cells infects the proportion of available uninfected cells. Increasing hf causes the proportion of available uninfected cells to become infected quicker and so exhausting the supply of cells to infect sooner (figure 3f). If this parameter is set to zero, then cell-to-cell infection cannot occur. Therefore, the rate of infection of the proportion of cells by free virions in the interstitial space (hy) must be sufficiently high to ensure that enough of the cells become infected to increase virus in the system. If hy is too low, then the virus exits the system via the interstitial space quicker than it can be produced and so decreases to zero. The proportion of infected cells, F, increases exponentially provided that cells first become infected from the interstitial space, and, second, the proportion of available cells for infection does not decrease too low.
As above, we calculate typical values for the cell-to-cell infection rate of the proportion of available epithelial cells for infection, given by the expression, hf lf FU, as shown in table 3. The overall infection rate given by the first two terms of dF/dt, (hy(Y2/(m2 + Y2))YU + hflfFU), is given in table 4. This shows that infection from the interstitial space contributes 50–60 per cent.
The parameter m represents a virus level below which there is a much reduced infectious ability, delaying the production of virus in the cell system (figure 7g). The maximum value of the 95% CIs of m for each of the high-dose animals is similar to one another. This is due to all of these datasets having few or no values below the detection threshold of 102 genomesml−1 (figure 6). Since virus levels in the high-dose animals are high, the reduced infectious term, Y2/(m2 + Y2), has less of an effect than it does for lower dose animals with lower virus levels. This is why higher doses have a larger confidence interval for m.
The common fit model shows that m is a necessary parameter, and hence justifies the inclusion of the term Y2/(m2 + Y2), having a lower 95% CI greater than zero. The value of m for the common fit model enables the correct dynamics to be displayed for all animals by only varying the initial condition. The parameter, m (genomesml−1), is estimated by the 95% CI for the common fit model: (2.11×10−1, 1.68×102).
The sensitivity to the constant lf apparent in figure 7h is considerable, but these graphs correspond to cell attachment and intracellular infection times of 24.24, 12.12, 6.06, 1.52, 0.76 and 0.38 hours for (1/8)lf, (1/4)lf, (1/2)lf, 2lf, 4lf and 8lf, respectively. These times are outside the calculated approximate time of between 2.5 and 3.5 hours.
When the amount of virus released per proportion of cells is increased, the infection rate of the proportion of cells also increases, consequently, exhausting the available cells for infection more quickly (figure 7i).
Decreasing the initial dose delays the production of virus, causing a viral load peak of the same magnitude later (figure 7j).
The parameter X(0) controls the initial quantity of virus in the system. A higher initial concentration of FMDV in the blood results in an earlier peak of virus in both the data and in the model. The model peak always has the same dynamics for any value of X(0), which can be seen clearly from the common fit graphs in figures 4 and and77d. From the 95% CI plots of X(0) from the individual fits, as shown in figure 5a, we can see that this fitted parameter is in good agreement with the initial dose given to each animal. In particular, it reflects that animals UQ10–UQ13 were given a low dose, animals UQ14–UQ17 a medium dose and animals UQ18–UQ21 a high dose. This is the only parameter that seems to reflect the fact that different animals were given different doses. From the 95% CI plots of X(0) from the common fit model, as shown in figure 5b, we also see the relationship with the initial dose. Although there is no clear difference between the low- and medium-dose animals, the timing of the peaks is similar.
From the graphs in figure 4, it appears that the model predicted peak values of X are identical for each animal. When the peak value is calculated for each animal, we see that all the peak values are equal at 1.79×108 FMDV genomesml−1. Unfortunately, it is not possible to write down an analytical solution for the peak value of X in terms of the given parameter set, since we must solve the system of differential equations numerically. However, in all the analyses performed in producing the best-fit common parameters, the peak value has always been the same across all animals.
In conclusion, although there is some variation in the dynamics of virus spread between animals, there is a general pattern that can be explained in terms of the initial dose of FMDV in the blood. The model predicts that, when the initial dose is varied, only the timing of the virus spreading is varied, with the dynamics remaining the same: the higher the initial dose the earlier the development of detectable viral load. This was shown by fitting the model to all 12 animal datasets only allowing the initial dose to vary between animals, as shown in figure 4.
The use of a model was important in this exercise providing estimates for missing data and providing a mechanistic explanation for the observed data. Furthermore, estimates for the parameters governing the within-host dynamics of FMDV were obtained. Some of the authors are currently involved in the ongoing experimental work to estimate the infectious and incubation periods of animals with FMDV, together with the implications for vaccines, an important component of this is the study of within-host dynamics.
Animal experimentation was approved by the Institute for Animal Health (IAH) ethical review board under the authority of a Home Office project licence in accordance to the Home Office Guidance on the Operation of the Animals (Scientific Procedures) Act of 1986 and associated guidelines.We would like to thank the BBRSC for funding this research.
In the model, the blood, X, and interstitial space, Y, are connected by a constant flow from X to Y with rate kxy. Similarly, there is a constant flow from Y to X with rate kyx, where r acts as a conversion factor for the concentration of FMDV between the different-sized spaces. Considering only the interactions between X and Y we have the following:
These equations have a set of steady states (i.e. X and Y remain constant when dx/dt and dy/dt equals zero) given by the line . It is easily shown that this set of steady states is stable. The equations have the following solution:
Therefore, as t→∞, and . Furthermore, since , we know that X and Y converge exponentially fast towards the stationary point (which is on the line kyxY=rkxyX) along the line Y=r(X(0)−X)+Y(0).
When we consider X and Y in the model, we observe that Y has extra interactions that increase and decrease the value of Y. However, the interactions of Y with X remain the same; therefore, when Y changes, X converges exponentially fast to the corresponding value on the line . This can be observed in figure 2, where the values of X and Y track one another.
From these observations of X and Y, we can conclude that it would be sufficient to model X and Y with only one variable. For example, we can drop the variable X from the model and achieve the same behaviour. The initial values of Y are then given by the original model initial values of X as .