Improving our forecasts for trachoma elimination: What else do we need to know?

The World Health Organization (WHO) has targeted trachoma for elimination as a public health concern by 2020. Mathematical modelling is used for a range of infectious diseases to assess the impact of different intervention strategies on the prevalence of infection or disease. Here we evaluate the performance of four different mechanistic mathematical models that could all realistically represent trachoma transmission. We fit the four different mechanistic models of trachoma transmission to cross-sectional age-specific Polymerase Chain Reaction (PCR) and Trachomatous inflammation, follicular (TF) prevalence data. We estimate 4 or 3 parameters within each model, including the duration of an individual’s infection and disease episode using Markov Chain Monte Carlo. We assess the performance of each models fit to the data by calculating the deviance information criterion. We then model the implementation of different interventions for each model structure to assess the feasibility of elimination of trachoma with different model structures. A model structure which allowed some re-infection in the disease state (Model 2) was statistically the most well performing model. All models struggled to fit to the very high prevalence of active disease in the youngest age group. Our simulations suggested that for Model 3, with annual antibiotic treatment and transmission reduction, the chance of reducing active disease prevalence to < 5% within 5 years was very low, while Model 2 and 4 could ensure that active disease prevalence was reduced within 5 years. Model 2 here fitted to the data best of the models evaluated. The appropriate level of susceptibility to re-infection was, however, challenging to identify given the amount and kind of data available. We demonstrate that the model structure assumed can lead to different end points following the implementation of the same interventions. Our findings are likely to extend beyond trachoma and should be considered when modelling other neglected tropical diseases.

The equations for Model 3 [1] are as follows: ∂I i ∂t ∂a i ∂t = a λ a S i (a, t) − µI i (a, t) − σI i (a, t) ∂ID i ∂t ∂a i ∂t = σI i (a, t) − µID i (a, t) − ω i ID i (a, t) + a (λ a P D i (a, t))Γ (7) ∂D i ∂t The equations for Model 4 are as follows: ∂I i ∂t ∂IO i ∂t ∂ID i ∂t ∂a i ∂t = η i I i (a, t) − µID i (a, t) − ω i ID i (a, t) + a (λ a D i (a, t))Γ (13) For all models evaluated the force of infection was structured in the manner, therefore λ is defined as: It has recently been suggested that the dynamics of trachoma may facilitate elimina-tion [2]. Here when we implement the different interventions, we assume that individuals are exposed to a linear a j α i ID i (a ,t) i N i (a ,t) and non-linear v 2 ( a i α i ID i (a ,t) i N i (a ,t) ) ψ + 1 force of infection (as previously modelled). For the non-linear term, at low prevalence it behaves similarly to a mass-action model , but at mid-range prevalence levels transmission will exhibit either positive or negative feedback depending on the value of v 2 [2]. Typically in this system positive feedback allows for two stable equilibria in a deterministic model, one with zero prevalence and one with a non-zero prevalence.
α i is the infectivity of an individual in compartment ID i . σ is the rate at which infected individuals I i progress to becoming infected and infectious, ω i is the rate at which ID i individuals recover from infection, ρ i is the rate at which D i individuals recover from active disease only. For Model 4 we assume that both IO and ID individuals contribute to the FOI so both are in the numerator for equation 15.
w(a, a ) is a mixing matrix which contains information on the rate of mixing between individuals of age group a and a' [3]. δ ( a, a ) is the Kronecker Delta. Kronecker Delta is a piecewise function of variables, if the two variables (in this instance age-group) are equal the value of the function is 1, otherwise it is 0. N a is the number of individuals of age a' in the population and indicates the degree of mixing assortativity, which can range between 0 (random age mixing) and 1 (fully assortative).

Rate of recovery
Rates of recovery and infectivity were modelled as previously described by Gambhir et al [4] The per individual rate of recovery ρ i from an individuals ith infection (measured as a rate per year). The rate of recovery is assumed to change as an exponential function of i that begins at a rate of ρ 1 (recovery from the first infection) and rises to a maximum rate ρ 1 00 where no greater rate of recovery can be achieved after this point [5].
We also model the per individual rate of recovery from active disease only ω i (measured as a rate per year) from infection i the rate of recovery is assumed to change as an exponential function of i that begins at a rate of ω 1 (recovery from the first infection) and rises to a maximum rate ω 100 . We assume the rate of change of the recovery rate per infection (φ) and per each disease episode (θ) were different for recovery from infection and active disease.
For the IO and PD states, we assume immunity follows the same functional form as above. We assumed the minimum duration in IO was the same as ID, and the maximum duration was half of the duration of time spent infectious estimated in [5] Values of these parameters are provided in Table S1.

Infectivity
We assumed the infectivity of an individual α i was proportional to their bacterial load.
In the model we reflect declines in bacterial load with repeated infection as reductions in the probability of transmission. We assumed that an individuals load (i.e probability of transmitting to others in the community) decreased with an increasing number of infections experienced. Therefore, infectivity was a function of the number of previous infections experienced by each individual, a trend that is in agreement with the data from trachoma endemic communities in which the bacterial load decreases with age [6][7][8]. We assumed a linear decline in infectivity for the 3 parameter models, and an exponential decline in The rate of decline in infectivity for the exponential decline in infectivity was: Model parameters and state variables

Model fitting
All parameter estimates in the 3 and 4 parameter model (β, α 1 , ρ 1 , ω 1 ) (Table S1)  For each parameter set we calculate the log-likelihood (LL), as previously done [4] of observing the age-specific data on prevalence of infection and disease, as well as the rates of recovery and infection load. We assumed a binomial likelihood for the prevalence data on infection and disease.
where p a is the model-predicted infection or disease prevalence for age group a; x a is the number of individuals that test positive for infection or disease and n a the total number of individuals in the data set for age group a.
To estimate the rate of recovery from an individuals first infection and disease episode as well as the rate of change in infectivity (as determined by infection load) we evaluated a Gaussian likelihood for each parameter separately.
Where V ar a is the age-specific fixed variance. Variance of the bacterial load and infection and disease durations for each age group was fixed due to the issue arising that some parameter sets that resulted in poorer fits to the data but minimised the variance resulted in better performing likelihoods. Therefore it was chosen to keep the variance constant across all parameter sets. The variance was fixed when the 3 and 4 parameter models were fit and for the iterations of the MCMC chain within each model. Variance was 128, 98, 32 days for each age group respectively for the duration of an infection and disease episode. Higher variance was assumed for the durations of infection and disease for the younger age group reflecting large uncertainty in previous estimates of these parameters [5]. For bacterial load a variance of 0.005 was used across all age groups. D a is the value in the data and M a is the model predicted value for a given set of parameters.
Therefore the 3 and 4 parameter models, the likelihood of each parameter set was eval-uated as the sum of each of the individual likelihood for each parameter.

Modelling interventions
In the analysis two interventions were modelled. Firstly we modelled the implementation of Facial cleanliness and Environmental improvements (F and E) as a reduction in the transmission rate parameter β. We considered 0, 10, 30 and 50% reductions in the transmission rate parameter over the intervention period.
We also performed annual rounds of mass drug administration. The total number of people effectively treated in each annual round was calculated as c * e, and this treatment was performed for all individuals who were in a PCR positive state. Treatment was assumed to have no impact on an individuals disease state. Therefore, members of the population that were treated successfully in the I i state were returned to the S i , thus we assumed that no improved immunity developed as a result of this infection. Individuals that were in the ID or IO states who cleared their infection following treatment progressed to the D state, hence they non-longer had a PCR detectable infection but were still TF positive. We assumed that immunity to infection did develop if an individual was treated in the ID or IO states. A schematic of this movement due to treatment is illustrated in Figure 1.

Results: Model fits and parameter estimates
Age specific fits to the bacterial load data are presented in Figure S1 and Table S4.
MCMC diagnostics to assess convergence of 2 MCMC chains using the Gelman-Rubin statistic (R c ) and the Effective Sample Size (ESS) of each posterior (Table S5). Ensuring    Figure S1: Estimates from the best performing models of the age-specific bacterial load data. Estimates of age-specific bacterial load data from statistically the best performing models for each structure evaluated. Data is shown in red,