An ODE model of yaws elimination in Lihir Island, Papua New Guinea

Introduction

Yaws is an infectious disease spread by skin to skin contact mostly amongst children (Marks, Solomon & Mabey, 2014). It is caused by bacteria Treponema pallidum ssp. pertenue, and begins at an abrasion or open wound which then develops into the “primary papule” or “mother yaw” (Perine et al., 1984). This initial stage is known as primary yaws and the lesion persists for three to six months (Marks et al., 2015a). A short latency period may occur after primary yaws if the primary papule naturally heals before secondary lesions develop (Perine et al., 1984). Secondary yaws begins with the appearance of other lesions anywhere on the body (Mitjà, Asiedu & Mabey, 2013). These lesions heal spontaneously resulting in a noninfectious latent period that, in some cases, may last the remaining lifetime of the person (Perine et al., 1984). During the latent period, previously infected individuals may relapse into secondary yaws up to 5 years after recovering from infection (Marks et al., 2015b). Up to 10% of yaws cases may develop into late yaws (Mitjà, Asiedu & Mabey, 2013), also known as tertiary yaws. Tertiary lesions tend to be very harmful with massive necrotic tissue destruction; yet they are noninfectious (Marks, Solomon & Mabey, 2014).

Yaws was the first disease to be targeted for eradication by the World Health Organization (Marks et al., 2015a). Before this initiative, 90 countries were reported as endemic, totaling about 50 million cases worldwide (Kazadi et al., 2014). The mass screening and treatment programmes led by WHO reduced the global prevalence by >95% between 1950 and 1964, but yaws has reemerged as a public health problem (Asiedu et al., 2008). In 2018, Papua New Guinea and the Solomon Islands reported over 10,000 suspected cases each (WHO, 2018b). According to the most recent report (WHO, 2020), 15 countries are considered as currently endemic for yaws; 87,877 suspected yaws cases were reported to WHO in 2020 from 11 countries from which 81,369 cases were in Papua New Guinea. Solomon Islands reported 13,694 cases in 2019. Figure 1 shows the most recent status of yaws endemicity.

A single oral dose of azithromycin was shown to be just as effective as the previous treatment of injectable penicillin (Mitjà et al., 2012). This initiated a new wave of interest in eradication of yaws (Marks et al., 2015a) and, in 2012, the WHO implemented the Morges Strategy to combat yaws transmission with the goal of eradicating the disease by 2020 (WHO, 2012). The Morges Strategy includes one or more rounds of total community treatment (TCT) where treatment is given to all members of the community and followed by the total targeted treatment(TTT) where treatment is administered to all actively infected individuals and their close contacts as a response to a local outbreak (WHO, 2012). This strategy continues to be the primary plan to eradicate yaws, but the timeline has been stretched to eradication by 2030 after the original 2020 goal was not met (Dyson et al., 2019).

Mathematical modelling is now a standard and indispensable tool for understanding disease dynamics and control (Anderson & May, 1992), yet there is a surprising lack of models of yaws transmission. Until 2012, the only math model of yaws considered the effect of the chicken pox virus on yaws (Gart & De Vries, 1966). More recently, the model from Mushayabasa et al. (2012) divides the population into rich and poor and preforms a theoretical analysis of disease-free and endemic equilibria. In Muench (2013), the authors fitted a simple catalytic model to age structured yaws data. All other models of yaws are stochastic and were designed to estimate various aspects of yaws eradication. In Fitzpatrick, Asiedu & Jannin (2014), the authors were concerned with the economic side of eradication and concluded that the eradication would not be expensive; yet there is still a large degree of uncertainty for the lack of available data (Dyson et al., 2019). In Dyson et al. (2017), the authors created a model to estimate the fraction of individuals that are missed during treatment in the Morges strategy. In modeling effectiveness of the Morges Strategy, Marks et al. (2017) investigated the probability of eradication. Fitzpatrick et al. (2018) created a linear regression model in order to predict the probability of case reporting in different previously affected countries based on different parameters. Mooring et al. (2019) builds off of the work of Marks et al. (2017), using the same compartment model and modeling the effects of different combinations of TTT and TCT. The most recent model comes from Holmes et al. (2020) in which the authors adapted the model from Dyson et al. (2017) to again simulate different combinations of TTT and TCT. It was found that different populations require different treatments, but in general, TCT was more effective in eradication.

The stochastic models such as recent Marks et al. (2017); Mooring et al. (2019); Holmes et al. (2020) are generally more suitable for the eradication end game than the deterministic compartmental models. Yet, the deterministic models are typically simple and easy to analyze, while still reasonably accurate and realistic. Given the lack of deterministic models of yaws transmission in general, our goal is to develop a deterministic model of yaws transmission and then use the model to compare the effectiveness of TTT and TCT strategies. We use the model to derive a formula for the basic reproduction number and to obtain simulated times needed for yaws elimination. Our model can be used as a quick estimate of the effectiveness of a particular treatment strategy.

Methods

We created a compartmental model shown in Fig. 2. Individuals are born as susceptible (S) at rate Λ. The susceptible individuals become exposed (E) after coming in contact with individuals having primary (Y₁) or secondary (Y₂) yaws; the transmission rate is β. After an incubation period lasting a time σ⁻¹, the exposed individual develops primary yaws and becomes infectious. The primary yaws lasts a time ${\lambda }_1^{-1}$, after which the individual may either develop secondary yaws (Y₂) with probability p_Y1Y2, or go into a first latency period (L₁) with probability p_Y1L1 = 1 − p_Y1Y2. This means that the rate of progression from Y₁ to L₁ is p_Y1L1λ₁ while the rate of progression from Y₁ to Y₂ is p_Y1Y2λ₁.

The average duration of the first latency period is ${\rho }_1^{-1}$; after that an individual develops secondary yaws. The average duration of secondary yaws is ${\lambda }_2^{-1}$. After this period, one could either develop rare, debilitating and very painful, but non-infectious, tertiary yaws (Y₃) with probability p_Y2Y3, or go into the second latent period with probability p_Y2L2 = 1 − p_Y2Y3. The average length of the second latency period is ${\rho }_2^{-1}$. Afterwards, individuals can relapse into secondary yaws with a probability p_L2Y2, or develop tertiary yaws with a probability p_L2Y3 = 1 − p_L2Y2.

Individuals can be treated and return to susceptible at a rate of τ_I for individuals in a compartment I ∈ {E, Y₁, Y₂, Y₃, L₁, L₂}. Treatment of each compartment depends on the elimination strategy and the specific values are discussed below.

Finally, all individuals are assumed to die at rate µ.

The parameters are summarized in Table 1. Most parameter values were estimated directly from the literature. The only two exceptions are the transmission rate β and the treatment rates τ. Details are shown in Appendix B.

The transmission rate β was obtained by fitting the endemic equilibrium to baseline data (prior the mass treatment, i.e., when all τ_I = 0) from Lihir Island in Papua New Guinea (Mitjà et al., 2015b).

To model TCT, we assume τ_I = 1/6 for all I, corresponding to treating the whole population every six months. To model TTT, we assume the best case scenario, i.e., τ_E = τ_Y1 = τ_Y2 = τ_Y3 = τ_L1 = 1/6 while τ_L2 = 0.1/6, i.e., we assume that the TTT strategy finds and treats only 10% of secondary latent cases every six months but otherwise finds and treat every other infected individual. This again corresponds to treating active yaws cases and all their closed contacts (that will be either exposed or at most the first latency period) once in six months. We assume that 90% of secondary latent cases are omitted in the treatment because the were infected independently many months or even years ago and are not close contacts to the currently acutely infected individuals. We adopted these assumptions since we can then demonstrate that even these high coverage, yaws will persist in the population for long time under TTT strategy. The protocol of Mitjà et al. (2015b) study also included a 2 year period of non-strategic treatment. For that period, we assumed τ_E = τ_Y1 = τ_Y2 = τ_Y3 = τ_L1 = 1/24 with coverage as in TTT.

We validated our model on data from the mass treatment trial in Lihir Island (Mitjà et al., 2018); see Fig. 3. The fitted curve follows general trends of the data. However, the real data for latent infections exhibits oscillation with peaks and dips every 6 months and our simple model cannot exhibit such oscillations.

We used the compartmental model from Fig. 2 to create a system of ordinary differential equations. We found disease-free and endemic equilibria. Using the next generation matrix method (van den Driessche & Watmough, 2002), we found the basic reproduction number. We performed the stability analysis of the disease-free equilibria based on methods from van den Driessche & Watmough (2002) and Castillo-Chavez et al. (2002).

We did simulations in MATLAB, the code is made available in supplementary material. We adhered to responsible coding practices as outlined in Lucas et al. (2020).

The global uncertainty and sensitivity analysis by the partial rank correlation coefficients, PRCC was based on Marino et al. (2008). We randomly selected 1000 parameter values within the ranges specified in Table 1. We used only those values that could fit to baseline data from Lihir Island (Mitjà et al., 2015b).

Results

We obtained an explicit formula for the basic reproduction number, R₀. As shown in Eq. (9), (1)${\displaystyle R_0=\left(\frac{\beta \sigma }{v_E{v}_{Y_1}}\right)\left(1+\left(\frac{{\lambda }_1{v}_{L_2}}{v_{L_1}}\right)\left(\frac{p_{Y_1{L}_1}{\rho }_1+p_{Y_1{Y}_2}{v}_{L_1}}{v_{L_2}{v}_{Y_2}-p_{Y_2{L}_2}{\lambda }_2{p}_{L_2{Y}_2}{\rho }_2}\right)\right)}$ where v_I denote the sum of all total rates out of the compartment I, i.e., (2)${\displaystyle v_E=\sigma +{\tau }_E+\mu }$ (3)${\displaystyle v_{Y_1}={\lambda }_1+{\tau }_{Y_1}+\mu }$ (4)${\displaystyle v_{Y_2}={\lambda }_2+{\tau }_{Y_2}+\mu }$ (5)${\displaystyle v_{Y_3}=\mu +{\tau }_{Y_3}}$ (6)${\displaystyle v_{L_1}={\rho }_1+{\tau }_{L_1}+\mu }$ (7)${\displaystyle v_{L_2}={\rho }_2+{\tau }_{L_2}+\mu .}$ We estimated that without any treatment, R₀ = 1.2548. The uncertainty analysis showed that in order to fit data from Lihir Island, R₀ is between 1.24 and 1.27; see Fig. 4. With TTT treatment, the values of R₀ ranged between 0.02 and 0.08. This shows that TTT is quite effective in prevention of the spreading of the epidemics.

We proved (Theorem 1 in Appendix A.2) that the disease-free equilibrium is globally asymptotically stable when the basic reproduction number R₀ < 1. We also showed (Lemma 1 in Appendix A.2) that R₀ is decreasing in the treatment rate. If the treatment rate is high enough, the basic reproduction number drops below 1 even for a conservative TTT strategy when only active cases of yaws gets treated; see Fig. 5. This means that the Morges strategy can eventually eliminate yaws.

To understand how long the Morges strategy needs to be applied, we simulated two rounds of initial TCT and followed by subsequent rounds of TTT. We performed a global uncertainty analysis where we varied parameters within the ranges specified in Table 1. Figure 6 demonstrates the results. Our model predicts that it would take about 14 to 16 years to achieve a thousandfold decrease in cases (i.e., less than 1 infected person in Lihir Island). The relatively high prevalence of latent cases in the population and the long latency period are the main culprits behind this long elimination time. The continuous application of TCT strategy every six months can achieve the same results in about 3.5 years; the improvement in speed is caused by the latent cases getting treated as well.

As illustrated in Fig. 7, the success or failure of TTT strategy significantly depends on how many latently infected individuals can be discovered and treated. The figure in fact shows expected elimination times for a whole family of strategies with TTT on one end (when the coverage of L₂ is low) and TCT on the other end (when the coverage of L₂ is 100%). It can take over 25 years to eliminate yaws if only 1% or less of latent cases are found; it would take about about 10 years if 20% is found and about 5 years if about 50% of cases is found.

The sensitivity analysis shows a strong influence of the relapse rate, ρ₂, and the spontaneous healing rate of the secondary yaws, λ₂, on the elimination time under the TTT regime; see Fig. 7. The higher the relapse rate and the lower the healing rate, the less it takes to eliminate yaws. This initially counter-intuitive result is caused by the fact that spontaneous healing increases the pool of latently infected individuals that can be missed by the TTT strategy. However, if infected individuals do not heal spontaneously, they can be discovered and treated. This again indicates that the latent individuals are the weakest point of the TTT strategy. A higher birth rate also reduces the time to elimination. This is mainly because a higher birth rate increases the influx of healthy individuals while the active yaws of young children are caught on time before progressing to latency. Naturally, a shorter incubation periods increases the time needed for the elimination as they increase the number of yaws cases. The effects of other parameters are relatively mild and not significant.

Finally, let us note that when R₀ > 1, the disease-free equilibrium is not stable and there exists an endemic equilibrium given explicitly in Eq. (10). We run numerical simulations for parameter values with ranges in Table 1 and the numerical solutions of the ODE model always converged to the endemic equilibrium. Moreover, motivated by Yang et al. (2017) and LaSalle (1976), we considered a Lyapunov function $\mathcal{L}={\sum }_C\left(C-C^{\ast }-C^{\ast }ln\left(\frac{C}{C^{\ast }}\right)\right)$, where the summation is taken over all compartments C ∈ {S, E, Y₁, Y₂, Y₃, L₁, L₂} and C^∗ is an endemic equilibrium value. It follows that $\mathcal{L}\ge 0$ and $\mathcal{L}=0$ iff C = C^∗ for all compartments. Also, we evaluated ${\mathcal{L}}^{\prime }={\sum }_C\left(1-\frac{C^{\ast }}{C}{C}^{\prime }\right)$ at 10⁵ randomly selected values of the compartments. We always saw that ${\mathcal{L}}^{\prime }<0$. Thus, we believe that the endemic equilibrium is globally stable whenever R > 1, although we do not have an analytical proof of this fact. However, as it has been shown in Fig. 4, even with the weaker TTT treatment, R₀ is significantly less than 1 and thus, for the purpose of the elimination (which is the main focus of this paper), the stability of the disease-free equilibrium is much more important.

Discussion

To model TTT strategy, we made a conservative assumption that not many latent cases are treated. We argue that this is a reasonable reflection of a reality in the eradication endgame. The latent cases represent reservoir of future infections Dyson et al. (2019). By treating a recently relapsed latent case with all its close contacts, TTT strategy prevents outbreaks. However, contact tracing does not identify many other latent cases in the population; they likely got infected independently many months or even years ago. Thus, TTT works quite slowly as an elimination strategy as it is equivalent to waiting for the latent cases to relapse instead of actively identifying and treating them while still asymptomatic.

Our model predicts very little variation of eradication times when using TCT strategy. This is natural as the whole population gets treated and most factors of yaws dynamics thus do not play any crucial role. The variability is much larger for the TTT regime which could potentially eliminate yaws in as little as 14 years but it may also take 16 years. The two key factors responsible for the large variation are the duration of the latent period (which is positively correlated with the elimination time) and the duration of the secondary yaws (which is negatively correlated with the elimination times). Gaining more knowledge about these two parameters would reduce the uncertainty of the model predictions.

Our model differs from previous models in two crucial aspects. First, we developed a deterministic ODE model, in contrast to recent stochastic models developed in Fitzpatrick, Asiedu & Jannin (2014); Marks et al. (2017); Dyson et al. (2017); Fitzpatrick et al. (2018); Holmes et al. (2020). While stochastic simulations can incorporate higher degrees of realism, there is a natural simplicity in the ODE models that allows for an easy estimation of the basic reproduction number. Even with different model parameters, we do not necessarily have to rerun the simulations to be able to predict the model outcomes. Our model can thus serve as a first and reasonably reliable estimate of what will happen under different elimination strategies.

Second, our model incorporates all of the known stages of yaws. All models of yaws should consider susceptible individuals, infectious stage(s) of yaws (possibly divided into primary and secondary yaws) and the asymptomatic latent yaws that can relapse. Similarly to Fitzpatrick, Asiedu & Jannin (2014), we also considered tertiary yaws; and as Mushayabasa et al. (2012), we included exposed individuals. Finally, as in Marks et al. (2017) we included a possibility of a latent period between the primary and secondary yaws.

There are several limitation of our model. Most of the limitations stem from the fact that our model is a simple deterministic ODE model in homogeneous population. It thus cannot capture the true eradication endgame when only very small, often just a single digit, number of individuals are infected. The model also cannot capture household dynamics as done in Dyson et al. (2017) or the population structured into hamlets as done in Mooring et al. (2019).

Unlike in stochastic simulations used in Marks et al. (2017); Dyson et al. (2017); Mooring et al. (2019); Holmes et al. (2020), we do not explicitly consider treatment coverage. An independent coverage is implicitly incorporated in our model—a rate τ_I = 1/6 can mean that the whole (100%) population is treated once every 6 months, as well as that the attempt to treat the whole population is made every m months but at each attempt, only p∗100% of the population is reached with m/p = 6. A systematic failure of the treatment could be included in the model by duplicating each compartment into “treatment adherent“ and “treatment non-adherent”. Setting the birth rate as (1 − p)Λ and pΛ, respectively, into the susceptible treatment adherent and treatment non-adherent, respectively, would then achieve a systematic failure of treatment for p∗100% of the population.

Economics plays a key role in the feasibility of yaws eradication. Our model should be extended by explicitly optimizing control strategies, i.e., the proper combination of TCT and TTT at the appropriate time intervals. The extensions need to take into the account that underdeveloped areas are more prone to transmission and are harder to screen for active infections.

Conclusions

Our paper is the first ODE compartment model specifically applied to yaws elimination. We investigated two strategies, the total community treatment (TCT) and the total targeted treatment (TTT). In agreement with previous models (Dyson et al., 2017; Marks et al., 2015c), we found that due to the high prevalence of latent infections, it is very hard to eliminates yaws by using TTT. Our model predicts that it would take about 15 years to reduce the prevalence thousandfold from the current levels. On the other hand, it would take only about 3.5 years if the whole community was treated once every six months. This is in a quantitative agreement with a recent detailed stochastic model (Holmes et al., 2020). We also note that due to the global stability of the disease-free equilibrium, and the fact that R₀ is significantly less than 1 even under TTT treatment, the initial levels of yaws in the population do not play a crucial role for the eradication.

In the light of above findings, we thus recommend using total community treatment as the primary yaws elimination strategy. This recommendation is further supported by the fact that (a) TCT provides additional benefits such as reduction in trachoma prevalence (Solomon et al., 2015), (b) the cost of TCT is not much larger that the cost of TTT (Fitzpatrick, Asiedu & Jannin, 2014), and (c) TTT requires active surveillance (Fitzpatrick et al., 2018), possibly further erasing the difference between the costs of these two approaches. As a note of caution, our model did not consider emergence of antibiotic resistant strains (Mitjà et al., 2018). It is a question whether a large scale application of TCT could eliminate yaws before the antibiotic resistance becomes a true obstacle.

Supplemental Information