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

View article
PeerJ

Main article text

 

Introduction

Methods

Results

Discussion

Conclusions

Supplemental Information

Code to draw a map

DOI: 10.7717/peerj.13018/supp-2

Additional Information and Declarations

Competing Interests

The authors declare there are no competing interests.

Author Contributions

Presley Kimball, Jacob Levenson, Amy Moore, Jan Rychtar and Dewey Taylor conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft.

Data Availability

The following information was supplied regarding data availability:

The MATLAB code for running the numerical ODE solutions is available in the Supplemental Files. The data used in this study is available in Table 1.

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.

Appendix A. Model analysis

 

The transmission dynamics described in Section yields the following system of ODEs dSdt=Λ+τEE+τY1Y1+τY2Y2+τY3Y3+τL1L1+τL2L2(βY1+Y2N+μ)S dEdt=βY1+Y2NS(σ+τE+μ)E dY1dt=σE(λ1+τY1+μ)Y1 dY2dt=pY1Y2λ1Y1+ρ1L1+pL2Y2ρ2L2(λ2+τY2+μ)Y2 dY3dt=pY2Y3λ2Y2+pL2Y3ρ2L2(μ+τY3)Y3 dL1dt=pY1L1λ1Y1(ρ1+τL1+μ)L1 dL2dt=pY2L2λ2Y2(ρ2+τL2+μ)L2.

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 N = S + E + Y1 + Y2 + Y3 + L1 + L2. The biologically feasible region consists of DR7+ such that NΛμ. Adding Eqs. (8)(14) yields dNdt=ΛμN and thus N(t)=Λμ(ΛμN(0))eμt.

Consequently, the region D is positively invariant and the model is epidemiologically and mathematically well-posed (Hethcote, 2000). We also see that limnN(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, E0=(S0,E0,Y01,Y02,Y03,L01,L02), is given by E0=(Λμ,0,0,0,0,0,0).

We calculate the basic reproduction number, R0, using the next generation method (van den Driessche & Watmough, 2002). Let the column vector I = (EY1Y2Y3L1L2)T represent the order of compartments with infection.

We define F and V as follows. The column vector F=(βY1+Y2NS,0,0,0,0,0)T represents new infections that are introduced into each compartment. The column vector V=((σ+τE+μ)E(λ1+τY1+μ)Y1σE(λ2+τY2+μ)Y2(pY1Y2λ1Y1+ρ1L1+pL2Y2ρ2L2)(μ+τY3)Y3(pY2Y3λ2Y2+pL2Y3ρ2L2)(ρ1+τL1+μ)L1pY1L1λ1Y1(ρ2+τL2+μ)L2pY2L2λ2Y2) 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., F=(0ββ000000000000000000000000000000000). Let V be the Jacobian of V at the disease-free equilibrium, i.e., V=(σ+τE+μ00000σλ1+τY1+μ00000pY1Y2λ1λ2+τY2+μ0ρ1pL2Y2ρ200pY2Y3λ2μ+τY30pL2Y3ρ20pY1L1λ100ρ1+τL1+μ000pY2L2λ200ρ2+τL2+μ). Let vI denote the sum of all total rates out of the compartment I, i.e., vE=σ+τE+μ vY1=λ1+τY1+μ vY2=λ2+τY2+μ vY3=μ+τY3 vL1=ρ1+τL1+μ vL2=ρ2+τL2+μ.

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 R0=(βσvEvY1)(1+(λ1vL2vL1)(pY1L1ρ1+pY1Y2vL1vL2vY2pY2L2λ2pL2Y2ρ2)).

Lemma 1

R0 ≥ 0. R0 is decreasing in τ. limτ→∞R0 = 0. Specifically, there is τ large enough such that R0 < 1.

Proof

1. Since vY2 = λ2 + τY2 + μ > pY2L2λ2 and vL2 = ρ2 + τL2 + μ > pL2Y2ρ2, all terms in Eq. (27) are non-negative. 2. We get R0τE=R0(σ+τE+μ)<0. Similarly, R0τY1=R0(λ1+τY1+μ)<0. Also, since τY2 appears only in the denominator of R0, R0τY2<0. It also follows that R0τL1<0 because τL1(pY1L1ρ1+pY1Y2(ρ1+τL1+μ)ρ1+τL1+μ)<0. Finally, R0τL2<0 because τL2(ρ2+τL2+μ(ρ2+τL2+μ)(λ2+τL2+μ)pY2L2λ2pL2Y2ρ2)<0.

More generally, it can also be shown that R0 decreases directly with respect to τ. Since dR0dτ=IR0τIτIτ and τI = cIτ, we have, by previous calculations, dR0dτ<0.

3. This follows directly from Eq. (27).□

Theorem 1

If R0 < 1, then the disease-free equilibrium is globally asymptotically stable. If R0 > 1, the disease-free equilibrium is unstable.

Proof

When R0 > 1, the disease-free equilibrium is unstable by van den Driessche & Watmough (2002).

When R0 < 1, the global stability follows from Castillo-Chavez et al. (2002). First, by Eq. (16), S0=Λμ, 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 = (EY1Y2Y3L1L2)T be the vector corresponding to infected compartments. The dynamics of I could thus be described by dIdt=(F − V)(S,I)=(FV)IˆG(S,I) where F and V are given by Eq. (19) matrix and Eq. (20) matrix and ˆG(S,I)=(β(Y1+Y2)(1SN),0,0,0,0,0)T. Note that F − V is an M-matrix (all the off-diagonal entries are non-negative) and all entries of ˆG(S,I) are non-negative since S ≤ N. Also, FV=DI(F − V)(S0,0). Thus, the assumption (H2) of Castillo-Chavez et al. (2002) is satisfied. Hence, the disease-free equilibrium (S00) is globally asymptotically stable when R0 < 1.□

Remark.

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

A.3 Endemic equilibrium

When solving for equilibria of the dynamics, i.e., the constant solutions of Eqs. (8)(14), we set the left-hand sides to zero and solve the following system of algebraic equations. 0=Λ+τY1Y1+τY2Y2+τEE+τY3Y3+τL1L1+τL2L2(βY1+Y2N+μ)S 0=βY1+Y2NS(σ+τE+μ)E 0=σE(λ1+τY1+μ)Y1 0=pY1Y2λ1Y1+ρ1L1+pL2Y2ρ2L2(λ2+τY2+μ)Y2 0=pY2Y3λ2Y2+pL2Y3ρ2L2(μ+τY3)Y3 0=pY1L1λ1Y1(ρ1+τL1+μ)L1 0=pY2L2λ2Y2(ρ2+τL2+μ)L2. Let N=S+E+Y1+Y2+Y3+L1+L2. By adding Eqs. (30)(36) and solving for N, we get N=Λμ. By Eq. (32), Y1=kY1E where kY1=σvY1. By Eqs. (35) and (39), 0=pY1L1λ1σvY1EvL1L1 and thus L1=kL1E where kL1=pY1L1λ1σvY1vL1. By Eq. (36), L2=pY2L2λ2Y2vL2. By Eqs. (33), (39), (41), and (42), 0=(pY1Y2λ1kY1+ρ1kL1)E+pL2Y2ρ2pY2L2λ2Y2vL2vY2Y2. Solving for Y2 yields Y2=kY2E where kY2=(pY1Y2λ1kY1+ρ1kL1)vL2vY2vL2pL2Y2ρ2pY2L2λ2 Note that, as in Lemma 1, kY2 > 0. By Eqs. (42) and (43), L2=kL2E where kL2=pY2L2λ2kY2vL2. By Eqs. (34), (43), and (44), Y3=kY3E where kY3=pY2Y3λ2kY2+pL2Y3ρ2kL2vY3. By Eqs. (31), (43), and (39), S=NvEβ(kY1+kY2).

Lemma 2

β(kY1+kY2)=vER0.

Proof

β(kY1+kY2)=βσvY1+βσλ1vL2(pY1Y2vL1+pY1L1ρ1)vY1vL1(vY2vL2pY2L2λ2pL2Y2ρ2) =(βσvY1)(1+(λ1vL2(pY1Y2vL1+pY1L1ρ1)vL1(vL2vY2pY2L2λ2pL2Y2ρ2))) =vE(βσvEvY1)(1+(λ1vL2vL1)(pY1Y2vL1+pY1L1ρ1vL2vY2pY2L2λ2pL2Y2ρ2)) =vER0.

Thus, by Lemma 2, Eqs. (47), (37), (39), (43), (45), (41), and (44), we have N=NR0+(1+kY1+kY2+kY3+kL1+kL2)E. The results of this section can be summarized in the following theorem.

Theorem 2

The endemic equilibrium E=(S,E,Y1,Y2,Y3,L1,L2) exists and is unique if R0 > 1. It is given by

S=ΛμR0, E=(Λμ)(11R0)(11+kY1+kY2+kY3+kL1+kL2) I=kIE,forI{Y1,Y2,Y3,L1,L2},

where kY1=σvY1 kL1=pY1L1λ1σvY1vL1 kY2=(pY1Y2λ1kY1+ρ1kL1)vL2vY2vL2pL2Y2ρ2pY2L2λ2 kL2=pY2L2λ2kY2vL2 kY3=pY2Y3λ2kY2+pL2Y3ρ2kL2vY3.

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 λ11=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 ρ11=1.5 months. All secondary yaws lesions subside in weeks to months (Mitjà, Asiedu & Mabey, 2013) and we will thus assume λ12=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 λ11=λ12=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 ρ12=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 pY2Y3 = 0.0001 and pL2Y2 = 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 pY1Y2 = 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, Y1 + Y2 + Y3, and L1 + L2 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.

Distribution of the parameter values used in the uncertainty and sensitivity analysis.

Figure 8: Distribution of the parameter values used in the uncertainty and sensitivity analysis.

The parameters were chosen with ranges as specified in Table 1 such that the model fits the data from Lihir Island (Mitjà et al., 2015b).

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.

References

 
1 Citation 1,281 Views 141 Downloads

Your institution may have Open Access funds available for qualifying authors. See if you qualify

Publish for free

Comment on Articles or Preprints and we'll waive your author fee
Learn more

Five new journals in Chemistry

Free to publish • Peer-reviewed • From PeerJ
Find out more