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

Yaws is a chronic infection that affects mainly the skin, bone and cartilage and spreads mostly between children. The new approval of a medication as treatment in 2012 has revived eradication efforts and now only few known localized foci of infection remain. The World Health Organization strategy mandates an initial round of total community treatment (TCT) with single-dose azithromycin followed either by further TCT or by total targeted treatment (TTT), an active case-finding and treatment of cases and their contacts. We develop the compartmental ODE model of yaws transmission and treatment for these scenarios. We solve for disease-free and endemic equilibria and also perform the stability analysis. We calibrate the model and validate its predictions on the data from Lihir Island in Papua New Guinea. We demonstrate that TTT strategy is efficient in preventing outbreaks but, due to the presence of asymptomatic latent cases, TTT will not eliminate yaws within a reasonable time frame. To achieve the 2030 eradication target, TCT should be applied instead.


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 Figure 1 World map of endemic history and prevalence of yaws. Data collected from WHO (2018b) and WHO (2020) and map was made with the aid of borders.m file (Greene 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 1 ) or secondary (Y 2 ) yaws; the transmission rate is β. After an incubation period lasting a time σ −1 , the exposed individual develops primary yaws and becomes infectious. The primary yaws lasts a time λ −1 1 , after which the individual may either develop secondary yaws (Y 2 ) with probability p Y 1 Y 2 , or go into a first latency period (L 1 ) with probability p Y 1 L 1 = 1 − p Y 1 Y 2 . This means that the rate of progression from Y 1 to L 1 is p Y 1 L 1 λ 1 while the rate of progression from Y 1 to Y 2 is p Y 1 Y 2 λ 1 . The arrows denote transitions between the compartments. The letters next to the arrows specify the per capita rates of the transitions. The red arrows denote a treatment. The black arrows show a natural disease progression (without any treatment), from susceptible (S) to exposed (E) and then to primary yaws (Y 1 ). After the primary yaws, a majority of cases goes through a short latent period (L 1 ) before progressing to secondary yaws (Y 2 ). A small portion develops secondary yaws directly. Most secondary yaws heals and the infection becomes latent (L 2 ). The latent cases can relapse to secondary yaws for the rest of their lives. Only a negligible number of individuals develop non-infectious tertiary yaws (Y 3 ).
Full-size DOI: 10.7717/peerj.13018/ fig-2 The average duration of the first latency period is ρ −1 1 ; after that an individual develops secondary yaws. The average duration of secondary yaws is λ −1 2 . After this period, one could either develop rare, debilitating and very painful, but non-infectious, tertiary yaws (Y 3 ) with probability p Y 2 Y 3 , or go into the second latent period with probability p Y 2 L 2 = 1−p Y 2 Y 3 . The average length of the second latency period is ρ −1 2 . Afterwards, individuals can relapse into secondary yaws with a probability p L 2 Y 2 , or develop tertiary yaws with a probability Individuals can be treated and return to susceptible at a rate of τ I for individuals in a compartment I ∈ {E,Y 1 ,Y 2 ,Y 3 ,L 1 ,L 2 }. 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 = τ Y 1 = τ Y 2 = τ Y 3 = τ L 1 = 1/6 while τ L 2 = 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  Asiedu & Mabey (2013) p Y 1 L 1 Probability of latency period after primary yaws Probability of immediate tertiary yaws infection after secondary yaws Asiedu & Mabey (2013) p Y 2 L 2 Probability of latency period after secondary yaws Probability of relapsing to secondary yaws during latent period after secondary yaws Asiedu & Mabey (2013) p L 2 Y 3 Probability of developing tertiary yaws during latent period See text 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  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 0 . As shown in Eq. (9), where v I denote the sum of all total rates out of the compartment I , i.e., We estimated that without any treatment, R 0 = 1.2548. The uncertainty analysis showed that in order to fit data from Lihir Island, R 0 is between 1.24 and 1.27; see Fig. 4. With TTT treatment, the values of R 0 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 0 < 1. We also showed (Lemma 1 in Appendix A.2) that R 0 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 2 is low) and TCT on the other end (when the coverage of L 2 is 100%). It can take over 25 years to eliminate yaws if only 1% or less of latent cases are  Table 1, τ 0,TCT ≈ 3.5 · 10 −4 and τ 0,TTT ≈ 1.68 · 10 −3 .
Full 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, ρ 2 , and the spontaneous healing rate of the secondary yaws, λ 2 , 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 We showed the dependence on the percentage of treated L 2 cases explicitly. The analysis of other parameters is done by partial rank correlation coefficients, PRCC (Marino et al., 2008).
Full-size DOI: 10.7717/peerj.13018/ fig-7 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 0 > 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 L = C C − C * − C * ln C C * , where the summation is taken over all compartments C ∈ {S,E,Y 1 ,Y 2 ,Y 3 ,L 1 ,L 2 } and C * is an endemic equilibrium value. It follows that L ≥ 0 and L = 0 iff C = C * for all compartments. Also, we evaluated L = C 1 − C * C C at 10 5 randomly selected values of the compartments. We always saw that L < 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 0 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) 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 0 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.

APPENDIX A. MODEL ANALYSIS
The transmission dynamics described in Section yields the following system of ODEs

A.1 Positivity and boundedness of solutions
All parameters in the model are non-negative and it can be shown that the solutions of the system Eqs. (8)- (14) are non-negative, given non-negative initial values. Indeed, let and thus Consequently, the region D is positively invariant and the model is epidemiologically and mathematically well-posed (Hethcote, 2000). We also see that lim n→∞ N (t ) = µ .
The system Eqs. (8)-(14) has two equilibria, the disease-free equilibrium and an endemic equilibrium as discussed below.

A.2 Disease-free equilibrium and the basic reproduction number
It follows from Eq. (17) general that the disease-free equilibrium, We calculate the basic reproduction number, R 0 , using the next generation method (van den Driessche & Watmough, 2002). Let the column vector I = (E,Y 1 ,Y 2 ,Y 3 ,L 1 ,L 2 ) T represent the order of compartments with infection.
We define F and V as follows. The column vector F = (β Y 1 +Y 2 N S,0,0,0,0,0) T represents new infections that are introduced into each compartment. The column vector represents difference between transfer out of the compartment and the transfer into that compartment that does not result from new infection. Let F be the Jacobian matrix of F at the disease-free equilibrium, i.e., Let V be the Jacobian of V at the disease-free equilibrium, i.e., Let v I denote the sum of all total rates out of the compartment I , i.e., Using symbolic computation capabilities of MATLAB, we can calculate V −1 , FV −1 and eigenvalues of FV −1 . The code is available in the supplementary material. Since only the first row of F is non-trivial, the same is true about FV −1 . Thus, FV −1 has only one non-zero eigenvalue given by More generally, it can also be shown that R0 decreases directly with respect to τ . Since and τ I = c I τ , we have, by previous calculations, dR 0 dτ < 0.3. This follows directly from Eq. (27).
Theorem 1 If R 0 < 1, then the disease-free equilibrium is globally asymptotically stable. If R 0 > 1, the disease-free equilibrium is unstable.
Proof When R 0 > 1, the disease-free equilibrium is unstable by van den Driessche & Watmough (2002). When R 0 < 1, the global stability follows from Castillo-Chavez et al. (2002). First, by Eq. (16), S 0 = µ , corresponding to the disease-free equilibrium, is globally asymptotically stable for Eq. (15). Thus, the assumption (H1) in Castillo-Chavez et al. (2002) is satisfied. Second, let I = (E,Y 1 ,Y 2 ,Y 3 ,L 1 ,L 2 ) T be the vector corresponding to infected compartments. The dynamics of I could thus be described by where F and V are given by Eq.  -Chavez et al. (2002) is satisfied. Hence, the disease-free equilibrium (S 0 ,0) is globally asymptotically stable when R 0 < 1.

Remark.
It follows from Castillo-Chavez et al. (2002) that there is C > 0 such that dI dt ≤ C(R 0 ) t . Consequently, the time needed the infections to drop below a predetermined level is proportional to (ln(R 0 )) −1 .

APPENDIX B. MODEL CALIBRATION
With the exception of the transmission rate β, the model parameters listed in Table 1 can be found directly in the literature. The birth rate in Papua New Guinea is 27.2 births per thousand per year (United Nations, 2019). The life expectancy of 65 years (World Bank, 2019). Should the model be used for other countries, we suggest to change these values appropriately since the uncertainty analysis suggest some sensitivity to these values as shown in Fig. 6.
The incubation period, σ −1 , after exposure to yaws lasts on average 21 days with a range from 9 to 90 days (Perine et al., 1984;WHO, 2018a). Primary lesions last for 3 to 6 months (Perine et al., 1984). We will assume λ −1 1 = 3 months since this allowed the best fit for our model; larger λ causes larger discrepancies between active and latent yaws cases. To estimate the length of the latent period after primary yaws, we note that secondary yaws occurs one to two months after the primary lesion heals (Marks et al., 2015a). We thus set ρ −1 1 = 1.5 months. All secondary yaws lesions subside in weeks to months (Mitjà, Asiedu & Mabey, 2013) and we will thus assume λ −1 2 = 3 months to be on par with the primary yaws. We note that we are primarily interested in the duration of the infectiousness. The estimated total duration of infectiousness for an untreated yaws patient, including relapses, is of the order of 12 -18 months (Perine et al., 1984). With λ −1 1 = λ −1 2 = 3 months, this would mean primary yaws, secondary yaws and two to four relapses into secondary yaws.
The second stage of latency ranges from zero to five years, and even up to ten years (Perine et al., 1984;Marks et al., 2015b). Thus, we assume ρ −1 2 = 30 months. Up to 10% of individuals develop tertiary yaws after five to ten years of untreated infection, but the condition is now extremely rare (Mitjà, Asiedu & Mabey, 2013). We thus set p Y 2 Y 3 = 0.0001 and p L 2 Y 2 = 0.9999. With these values, our model estimates the prevalence of the tertiary yaws in the endemic equilibrium (of untreated population) under 0.1% of the total population and slightly under 5% of the number of secondary yaws cases.
About 9-15% of primary yaws cases progress into secondary yaws with the primary lesion still present (Mitjà, Asiedu & Mabey, 2013;Marks et al., 2015b). We thus assumed p Y 1 Y 2 = 0.12 for the probability that individual progresses directly from the primary yaws to secondary yaws without any noticeable latent period.
To estimate the value of the transmission rate β, we fitted our model to the baseline data from Lihir Island (Mitjà et al., 2015b). The study indicates that, before the mass treatment trial, 2.4% of the population had active yaws and 18.9% were in the latent stage. The endemic equilibrium of our model given by Theorem 2 can be expressed as a function of β. We used MATLAB's optimization toolbox to numerically find the value of β so that the endemic equilibrium distribution of S + E, Y 1 + Y 2 + Y 3 , and L 1 + L 2 fits best the empirical values 0.797, 0.024, and 0.189. We got β = 0.016581. With this value, the data from Lihir Island were recovered with an error less than 0.005; most of the error was caused by underestimating active yaws cases.
We note that we ran a number of different scenarios and by allowing a non-zero probability for the individuals with primary or secondary yaws infections to completely recover from the infections and become susceptible, we were able to match data from Lihir Island with the error less than 10 −6 (and very similar model outcomes, in particular still 3.5 years to elimination by using TCT and about 20 years to elimination by using TTT).

APPENDIX C. GLOBAL UNCERTAINTY AND SENSITIVITY ANALYSIS
We performed the global uncertainty and sensitivity analysis by the partial rank correlation coefficients (Marino et al., 2008). For every parameter with the exception the treatment rate τ I (that does not have any influence at the base prior the beginning of the treatment) and the unknown transmission rate β, we randomly assigned a value from the uniform distribution within the range specified in Table 1. We then tried to fit the transmission rate β to the baseline data while keeping the treatment rate 0. We used the Matlab optimization toolbox as described in the previous section. We then used 1000 sets of the values that could fit to baseline data from Lihir Island (Mitjà et al., 2015b) within an error of 0.005 or less. The actual distribution of the parameter values used in the uncertainty analysis is shown in Fig. 8.
For each of these parameter values, we calculated the basic reproduction number, and the time needed for yaws elimination under the TCT and TTT protocols. The resulting histograms are shown in the main body of the manuscript.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
Presley Kimball, Jacob Levenson and Amy Moore were supported by the VCU REU program in mathematics funded by the National Security Agency grant number H98230-20-1-0011 and by the National Science Foundation grant number DMS1950015 awarded to Dewey Taylor. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.  Table 1 such that the model fits the data from Lihir Island (Mitjà et al., 2015b).