Stochastic SIS Model of Ocular Chlamydial Infection

Previously, we constructed a simple stochastic SIS model (Susceptible, Infected, Susceptible) of ocular chlamydial infection in a core group of children [20] (link),[21] (link). In this report, the previous model was modified to include infection from outside of the community, and a transmission term that can vary between communities in the same region as a normally distributed random effect (to account for the known variability of communities). Treatment strategies that allowed graduation of communities when observed infection fell below a certain threshold, as detected by a POC test, were incorporated. Specifically, we constructed a Markov model by letting denote the probability that there are i infected individuals in the population at time t (where i varies from 0 to N). At scheduled treatments, we assumed that each infected individual had an 80% chance of being treated (the WHO-recommended coverage rate), and if treated would revert to being susceptible. Between periodic mass treatments, the model is a standard continuous time Markov process. We assumed equilibrium at baseline, that infected individuals recover naturally at rate γ, that uninfected individuals become infected at rate β I/N from sources within the community (with β = R₀·γ), and at a rate of ν from outside the community (with ν decreasing to zero once wide-spread programs have begun), leading to the following set of N+1 Kolmogorov-forward equations:
For clarity, we expressed β in terms of the basic reproductive number, R₀ (where β = R₀·γ). Note that R₀ is defined as the mean number of secondary infectious cases caused by a single infectious case in an otherwise completely susceptible community [22] . At the time of the scheduled periodic mass treatments, we assume that each infected individual has a probability c of being treated (the effective coverage), with the number of infections post treatment being drawn from the corresponding binomial distribution.
Parameters for this stochastic model were fitted to baseline and 6-month data for each country using maximum likelihood estimation. We initiated simulations at the average prevalence for that region, and simulated the Kolmogorov-forward equations for 40 years to allow the distribution of prevalence to approximate the pre-treatment distribution at time point zero. We also initiated the model at the observed 2-month prevalence and simulated the equations for 4 months to estimate the expected distribution of prevalence at 6-months. The total log-likelihood was the sum of the baseline and the 6 month log-likelihoods for each of the communities in the area. Note that any event that occurred between baseline and 2-months (such as treatment, or mass re-infection from travel) would not bias these results [8] (link). Based on these Kolmogorov equations, the values of the parameters R₀, standard deviation of R₀ (thus treating R₀ as a random effect), γ, and ν that maximized the probability of obtaining the observed baseline and 6-month data for that country (i.e. the likelihood) were determined using an iterative, hill-climbing algorithm. Numerical optimizations were repeated a minimum of 4 times from random starting points (because of the possibility any single run could converge to a local, rather than the global, maximum); each iteration converged to the same value. Furthermore, a grid search did not reveal any greater maxima.
The variance of these estimates was assessed by inverting the Hessian of the log-likelihood evaluated at the maximum likelihood estimate (although note that the 95% confidence interval could not include ν = 0, because in each country, a community went from 0 infections at 2 months to >0 infections at 6 months). Coverage was assumed to be 80%, and the average population size was set at the mean of empirical results from the surveyed communities in that region (Table 1). For sensitivity analyses, we ran 1000 simulations with the fitted parameters under different scenarios. We varied the sensitivity and specificity of the POC test, as well as the threshold for declaring graduation, keeping R₀, standard deviation of R₀, γ, and ν at the optima found for that region. If not being varied, the sensitivity and specificity of a POC diagnostic test were set at 70% and 99% respectively, and the prevalence threshold for graduating communities of ≤5%. All analyses were carried out in Mathematica 5.2.