A stochastic model of antibiotic misuse, economy, and drug resistance: relating mutant extinction probability to socioeconomic and epidemiological factors

ABSTRACT

Controlling antibiotic drug resistance requires understanding extinction and emergence mechanisms of emerging bacteria. Selective pressure from prolonged antibiotic misuse may cause high-level antimicrobial resistance. Self-medication and other socioeconomic factors reduce antibiotic use, accelerating the emergence and extinction of resistant pathogens through stochastic fluctuation. This continuous antibiotic self-medication exposes individuals and communities to antibiotic resistance, especially in low- and lower-middle-income countries, according to current literature. We developed a stochastic drug-resistance model that integrates socio-economic growth and antibiotic use to study extinction and strain establishment in this paper. We analytically derived the extinction threshold using the multi-type branching process and obtained pathogen extinction conditions that match numerical experiments. The model's sensitivity analysis identifies extinction dynamics' key parameters. Our results show that higher income, awareness, and lower antibiotic use may increase the chance of extinction by reducing antibiotic misuse, along with strain transmission potential. These analyses may help public health policymakers combat drug resistance.

KEYWORDS:

• Drug resistance
• Stochastic model
• Socioeconomi growth

1. Introduction

The ever-increasing spread of drug resistance is a serious threat to global health, and it has received a lot of attention in the microbiological and medical communities [1]. The evolution of sensitive strain through mutation and plasmid transfer accelerates the drug-resistance process [2–4]. While the development of mutant variants is a random process that occurs at the within-host level, plasmids are transmitted from one bacterium to another (even from another species) at the population level, primarily through conjugation [5].

In order to establish themselves within the population, phenotypically different variations of mutant strains emerge in a small number of hosts and eventually escape from stochastic extinction and transmission in the population [6]. Endless mutant strains demonstrating various levels of resistance are thus able to rigidly fix themselves in the population [7]. For example, resistance occurs in Mycobacterium tuberculosis usually through point mutations, whereas many other bacterial pathogens can acquire resistance by horizontal transfer [8]. Mutation rates also vary widely – for example, the mutation rate of influenza is very high compared to any other pathogen, whereas it can be several orders of magnitude lower in the case of malaria than that of influenza [9; 10]. In contrast, vancomycin-resistant S. aureus or the spread of the gene NDM-1 develops resistance through the transfer of genetic material between species via horizontal gene transfer [11; 12]. Often plasmids may contain several resistance genes so that one drug can select for multi-resistance [13]. The rapid spread of plasmids through bacterial populations enables them to escape negative selection by moving to other hosts [14]. The dynamics of resistance under plasmid-borne and mutational mechanisms might therefore differ, but these differences are not yet largely explored.

On the other end of the spectrum, there is a strong association between the economic growth of nations and the current global drug-resistance scenario [15–17]. The ever-increasing antibiotic consumption and development of resistance are observed in many countries worldwide (Figures E1 and E2 in Appendix). A report by the World Bank Group in 2016 entitled Drug Resistant Infections: A threat to Our Economic Future indicates that resistance has the potential to cause similar or more economic damage than the 2008 financial crisis [18]. According to an expected future projection in this paper, medication resistance can lead low-income nations to lose more than 5% of their GDP and lock approximately 28 million people in poverty by the year 2050. The increase in the current global burden of resistance was mainly due to high consumption in LMIC through self-medication [19]. According to the World Health Organization (WHO), nearly 80% of the antibiotics used in LMICs are used in community settings, with approximately 20–50% of the antibiotics used inappropriately, while more than two-thirds of antibiotics available in the pharmaceutical sector are used for self-medication [18–20]. A new study shows self-medication incidence data in WHO Southeast Asian Region nations [21] (See Figure E2 in the Appendix). However, unregulated and improper antibiotic use has an impact on the building of resistance in the population by altering the pathogen’s appearance and extinction process.

In order to determine the probability of an outbreak or the disease extinction for such infections, a Continuous-Time Markov Chain (CTMC) framework is more efficient than a deterministic version of the model and its metrics $R_0$ [40–44]. The theory of multi-type branching process is used to determine the stochastic threshold of extinction [42; 45; 46; 47]. Many stochastic models have been developed to estimate the probability of extinction of infections and organisms applying branching processes in the field of ecological populations, cellular processes, genetics, and epidemics on networks [48–50]. However, to our knowledge, there is no modelling work emphasizing discussions on the stochastic extinction of drug resistance that incorporates the socioeconomic factors involving self-medication in the development of drug resistance in the community.

In our previous work, we developed a model to understand the dynamics of drug-resistance phenomena involving socio-economic growth in population [15]. The earlier model was built upon several simplistic assumptions such as linear economic growth over time, and the linear effect of income on the inappropriate use of antibiotics through self-medication. In addition, no mutation or plasmid transfer was assumed in the earlier model. We relaxed such assumptions and enhanced the model in the current framework. We establish the current framework as a CTMC model, which describes the random behaviour that appears in the transmission of resistant strain in the population.

We perform both analytical derivation and numerical computation of the extinction threshold, and they exhibit the right match. We also perform a sensitivity analysis to identify the parameters that have the most impact on the extinction of the strains. Our model demonstrates that the transmission potential of respective strains plays a major role in extinction, but at the same time, higher income, awareness, and lower antibiotic use may increase the chance of establishment and extinction significantly. These analyses are beneficial to health policy-makers and may quantify some parameters that are important to controlling the situation.

The organization of paper is as follows: In Section 2, we describe the integrated framework of antibiotic usage, socioeconomic growth, and transmission of drug-resistant strain in the population, Section 3 describes basic model properties, Section 4 illustrates the stochastic implementation of the dynamic model and computes the probability of extinction, Section 5 describes numerical computation and results. Finally, we discuss the limitations and scope of the research in this area.

2. Model

In this paper, we extend our previous deterministic drug-resistance model to a stochastic framework where emergence and selection of resistant strain, the spread of infection, and volume of drug use are highly stochastic in nature, although the skeleton of the stochastic model is based on an amalgamation of socio-economic growth, ecology of infectious disease, and antibiotic (mis)use by individuals including mutation and plasmid transfer [15]. See Figure 1 for different components and interactions among them. Here, we briefly describe different components of the model:

Figure 1. Schematic of the integrated model system combining (a) infection process component (b) Economic production component, and (c) Antibiotic overuse components. Arrows indicate feedback from one component on the other generating the self-reinforcing dynamics of the system.

2.1. Economic growth component

In the economic growth component, we used the Constant Elasticity of Substitution production function, introduced by Solow [51], that describes that the production technology has a constant percentage change in capital due to marginal change in technical substitution [51,52]. The effective output per unit labour (i.e. capital) $h(t)$ is given by the following equation:

(1) $g(t)={\left[r_{h}h{(t)}^{\frac{\varsigma -1}{\varsigma }}+r_l\right]}^{\frac{\varsigma }{\varsigma -1}}$(1)

where $g(t)$ is the production function, $r_h$ and $r_l$ are shared in the output from capital and labour respectively with $r_h+r_l=1$. $\varsigma$ is the elasticity of substitution between capital and technology-augmented labour, determining the homogeneity of the production function. Leontief, linear and Cobb Douglas functions are special cases of the CES production function [51]. Considering depreciation, the rate of change of capital becomes:

(2) $\frac{dh}{dt}=g(t)-{\delta }_{1}h(t)$(2)

where $h(t)$ is the capital (income) at time $t$ scaled to the labour supply and ${\delta }_1$ is the rate of capital depreciation.

2.2. Disease transmission component

We extend the classic Susceptible-Infected-Susceptible (SIS) framework to model disease transmission due to both strains – sensitive and resistant. We assume drug-resistant strain is already present and circulating in the population, and individuals are either colonized by sensitive strains or by drug-resistant strains and not both. The population is a closed population and divided into four compartments: Susceptible not using antibiotics $x(t)$, susceptible using antibiotics $S(t)$, infected by sensitive strain $y(t)$, and infected by resistant strain $z(t)$. We assume that a fraction of individuals colonized with sensitive strains will select pre-existing antibiotic-resistant mutants, and the resistant population increases proportionally to antibiotic usage (in unit Daily Defined Dosages (DDDs)). In contrast, plasmid transfer occurs through contacts between individuals colonized with a sensitive strain and colonized with a resistant strain at the rate $\xi$. The parameter $\mu$ defines natural mortality and $d$ is death due to the severity of drug resistance. The rate of mutation is given by $\sigma$. $1/{\gamma }_1$ is the average duration of colonization of sensitive strain and $1/{\gamma }_2$ is the average duration of colonization of resistant strain. $1/{\delta }_2$ is the average duration of antibiotic consumption. $p$ is the fraction of the recovered individuals who are still taking antibiotics. Figure 1 exhibits the schematic of the interaction among disease compartments. The following coupled differential equations govern the dynamics:

$\frac{dS}{dt}=a_{h}x-\beta { }^{\mathrm{\prime }}Sz+p{\gamma }_{2}z-{\delta }_{2}S+(1-\sigma )a_{h}y-\mu S$

$\frac{dy}{dt}=\beta xy-{\gamma }_{1}y-a_{h}y-\xi yz-\mu y$

(3) $\frac{dz}{dt}=\beta { }^{\mathrm{\prime }}\left(x+S\right)z-{\gamma }_{2}z+\sigma a_{h}y+\zeta yz-\mu z-dz$(3)

where $x(t)+S(t)+y(t)+z(t)=1$. We can exclude the differential equation of $x(t)$ from system (3), because population is a closed population.

2.3. Antibiotic misuse

Irrational use of antibiotics has become a frequent practice in developing countries. Self-medication and over-the-counter (OTC) drug use is widespread in low and lower-middle-income countries (LMICs) [53–57]. Countries in the WHO Southeast Asian region, like India, Bangladesh, Nepal, Bhutan, Indonesia, Korea, Sri Lanka, and Thailand, have high self-medication [21; 56] (see Figure E1 and E2 in Appendix). However, there is no empirical evidence against this complex nonlinear association between self-medication practice and the economic growth of the population. We consider an exponential function to represent the reciprocal relation between socioeconomic status and self-medication in a population compared to the linear relation assumed in our previous paper [15; 58]:

(4) $a_h(h)=(\hat{a}-\tilde{a})e^{-h{h}_o}+\tilde{a}$(4)

To test our hypothesis, we illustrate and compare the complex outcomes as emergent properties of parsimonious integrated models. In the exponential framework, antibiotic (mis)use is high, say $\hat{a}$ when income (capital) is low, and approach a minimum level $\tilde{a}$ whenever income is high.

Figure 2 exhibits the function $a_h$ schematic. The parameter $h_o$ is a crucial parameter that indicates the sensitivity of public’s education and awareness level to antibiotic utilization [15; 59]. The area under the curve determines antibiotic consumption. Figure 2 explains that the lower the value of $h_o$, the higher the self-medication, and hence higher the antibiotic consumption.

Figure 2. The schematic of the antibiotics consumption as function of income. The shaded region under the curve defines the volume of antibiotic use corresponding to the different values of $h_0$. From upper to lower curves the values of $h_o$ are taken as {0.5, 1, 2, 3}.

2.4. Integrated model system

We consider the share parameters $r_h$ and $r_l$ in the growth model (1) to be functions of strain prevalence because the most immediate impact of health on economic growth occurs through the cumulative effects of all infections [19; 20; 60; 61; 62]:

(5) $r_h(y,z)=r_h^c(1-y)(1-z)$(5)

(6) $r_l(y,z)=r_l^c(1-y)(1-z),$(6)

where parameters $r_h^c>0$ and, $r_l^c>0$ determine the growth of capital in the absence of any infection in the population and ($r_h(y,z)+r_l(y,z)\le 1)$. Combining (1) - (6), we have the following coupled differential equations representing the integrated system as follows:

$\frac{dS}{dt}=a_{h}x-\beta { }^{\mathrm{\prime }}Sz+p{\gamma }_{2}z-{\delta }_{2}S+(1-\sigma )a_{h}y-\mu S$

$\frac{dy}{dt}=\beta xy-{\gamma }_{1}y-a_{h}y-\xi yz-\mu y$

$\frac{dz}{dt}=\beta { }^{\mathrm{\prime }}(x+S)z-{\gamma }_{2}z+\sigma a_{h}y+\xi yz-\mu z-dz$

(7) $\frac{dh}{dt}=[r_h{h}^{\frac{\varsigma -1}{\varsigma }}+r_l{]}^{\frac{\varsigma }{\varsigma -1}}-{\delta }_{1}h$(7)

where $a_h(h)=(\hat{a}-\tilde{a})e^{-h{h}_o}+\tilde{a}$, and $x(t)+S(t)+y(t)+z(t)=1$.

3. Basic model properties

The initial conditions: $(x(0),S(0),y(0),z(0),h(0){)}^T\in$ $\left\{{(x,S,y,z,h)}^T\in R_{+}^5\right\}$ and the parameters are real positive. In order to show that our model is well-posed and biologically meaningful, we will prove that the solution set of the given model is positively invariant [63].

Lemma 3.1

Each solution $(S(t),y(t),z(t),h(t))$ of the model with the non-negative initial condition is non-negative for all $t>0$.

Proof.

Let $\mathcal{V}=(S,y,z,h{)}^T$

We can rewrite the model as

$\dot{\mathcal{V}}=f(\mathcal{V})$

where

$f(\mathcal{V})=\left[\begin{array}{c}{f}_1(\mathcal{V})\\ f_2(\mathcal{V})\\ f_3(\mathcal{V})\\ f_4(\mathcal{V})\end{array}\right]$

$=\left[\begin{array}{c}{a}_{h}x-\beta { }^{\mathrm{\prime }}Sz+p{\gamma }_{2}z-{\delta }_{2}S+(1-\sigma )a_{h}y-\mu S\\ \beta xy-{\gamma }_{1}y-a_{h}y-\xi yz-\mu y\\ \beta { }^{\mathrm{\prime }}(x+S)z-{\gamma }_{2}z+\sigma a_{h}y+\xi yz-\mu z-dz\\ {\left[r_h{h}^{\frac{\varsigma -1}{\varsigma }}+r_l\right]}^{\frac{\varsigma }{\varsigma -1}}-{\delta }_{1}h\end{array}\right]$

Now,

$\frac{dS}{dt}{|}_{S=0}=a_{h}x+p{\gamma }_{2}z+(1-\sigma )a_{h}y\ge 0$

$\frac{dy}{dt}{|}_{y=0}=0$

$\frac{dz}{dt}{|}_{z=0}=\sigma a_{h}y\ge 0$

$\frac{dh}{dt}{|}_{h=0}=r_l^{\frac{\varsigma }{\varsigma -1}}$

Following the invariant property of a dynamical system (see [63]), we see that $R_{+}^4$ is positively invariant set for the system given by equations (7)(7) $\frac{dh}{dt}=[r_h{h}^{\frac{\varsigma -1}{\varsigma }}+r_l{]}^{\frac{\varsigma }{\varsigma -1}}-{\delta }_{1}h$(7) . This implies that our problem is mathematically well-posed.

3.1. Basic reproduction ratio

The disease-free equilibrium (DEF) is given by: $\left(x^{\ast },S^{\ast },y^{\ast },z^{\ast },h^{\ast }\right)=\left(\frac{\left(\mu +{\delta }_2\right)N}{{\delta }_2+\mu +a_h},\frac{N{a}_h}{\mu +{\delta }_2+a_h},0,0,{\left(\frac{r_l}{{\delta }_1^{\frac{\varsigma -1}{\varsigma }}-r_h}\right)}^{\frac{\varsigma }{\varsigma -1}}\right)$

We use the next-generation matrix method to find out the basic reproduction number ${\mathcal{R}}_o$ [64–67]. We have two compartments $y$ and $z$ where new infections appear. So, calculating the next-generation matrix from a model (7), we have $K=F{V}^{-1}=\left[\begin{array}{cc}{\mathcal{R}}_{oy}A& 0\\ {\mathcal{R}}_{oz}B& {\mathcal{R}}_{oz}\end{array}\right]$

The spectral radius of $F{V}^{-1}$ i.e, $\rho (K)$ is the basic reproduction number which is given by:

(8) ${\mathcal{R}}_o=max\{{\mathcal{R}}_{oy}A,{\mathcal{R}}_{oz}\}$(8)

where ${\mathcal{R}}_{oy}=\frac{\beta }{\mu +{\gamma }_1}$, ${\mathcal{R}}_{oz}=\frac{\beta { }^{\mathrm{\prime }}}{\mu +p{\gamma }_2+d}$, $A=\frac{(\mu +{\gamma }_1)(\mu +{\delta }_2)N}{(\mu +{\gamma }_1+\sigma a_h)(\mu +{\delta }_2+a_h)}$ and $B=\frac{\sigma a_h}{\mu +{\gamma }_1+\sigma a_h}$ are transmission potentials of sensitive and resistant strains, respectively.

(1) If ${\mathcal{R}}_o\le 1$, then DFE is globally asymptotically stable.

(2) If ${\mathcal{R}}_o>1$, then there exists a unique endemic equilibrium which is locally asymptotically stable.

See appendix A1 for more information.

4. Stochastic implementation of integrated model

The basic reproduction number articulated above can be used to determine, on average, whether the appearance of a small number of infectious host pathogens will result in a major outbreak or disease will eventually die out. However, stochastic extinction may occur even if ${\mathcal{R}}_o>1$, which is called stochastic fade-out of an epidemic [44,46]. It is necessary to consider a stochastic model in order to determine the probability of disease extinction or predict a major outbreak that occurs by a pathogen. In this section, we implement the same deterministic model (7) to a continuous-time Markov chain (CTMC) model with a finite number of population and considered the multi-type Galton-Watson branching process (GWbp) to calculate the probability generating functions (pgf) from the CTMC model. For a general idea about the GWbp, please see the references [46,68,69]. It may be worthwhile to mention here that the conditions $x(t)+S(t)+y(t)+z(t)=1$ does not hold further for the stochastic model.

4.1. CTMC model and derivation of extinction probability

We implement $x(t),S(t),y(t)$, $z(t)$ and $h(t)$ as discrete-valued random variables at time $t$. The dynamics of the CTMC model can be interpreted in terms of the infinitesimal transition probabilities for stochastic processes using the vector $\overrightarrow{X}\in R_{+}^5$ given as

$\overrightarrow{X}(t)=(x(t),S(t),y(t),z(t),h(t)),t\in [0,\mathrm{\infty }).$

The transition probabilities and their rates are given in Table 1.

Table 1. State transitions and rates.

Description	State transition a$\to$b	Rate $P$(a,b)
Birth of	$(x+1,S,y,z,h)$	$\mu$ N
Birth of	$(x,S+1,y,z,h)$	$a_{h}x$
Birth of	$(x-1,S,y+1,z,h)$	$\beta xy$
Birth of	$(x,S-1,y,z+1,h)$	$\beta { }^{\mathrm{\prime }}Sz$
Recovery of	$(x+1,S,y-1,z,h)$	${\gamma }_{1}y$
Recovery of	$(x,S+1,y,z-1,h)$	$p{\gamma }_{2}z$
Birth of by mutation	$(x,S,y-1,z+1,h)$	$\sigma a_{h}y$
Birth of by plasmid transfer	$(x,S,y-1,z+1,h)$	$\xi yz$
Birth of by susceptibles	$(x-1,S,y,z+1,h)$	$\beta { }^{\mathrm{\prime }}xz$
Maximum infected treatable	$(x+1,S,y,z-1,h)$	$(1-p){\gamma }_{2}z$
Birth of who does not acquire resistance by mutation	$(x,S+1,y-1,z,h)$	$(1-\sigma )a_{h}y$
Duration of antibiotic use	$(x+1,S-1,y,z,h)$	${\delta }_{2}S$
Natural death of	$(x-1,S,y,z,h)$	$\mu x$
Natural death of	$(x,S-1,y,z,h)$	$\mu S$
Natural death of	$(x,S,y-1,z,h)$	$\mu y$
Natural death of	$(x,S,y,z-1,h)$	$\mu z$
Death of due to resistance	$(x,S,y,z-1,h)$	$dz$
Growth of	$(x,S,y,z,h+1)$	$[r_h{h}^{\frac{\varsigma -1}{\varsigma }}+r_l{]}^{\frac{\varsigma }{\varsigma -1}}$
Discounting of	$(x,S,y,z,h-1)$	${\delta }_{1}h$

Using the model equations in (7), we compute the probability generating function (pgf) for offspring from sensitive strain $y$, given $y(0)=1$ and $z(0)=0$ is

$f_1(u_1,u_2)=\frac{({\gamma }_1+\mu )+(\xi +\sigma a_h)u_2+\beta u_1^2}{\beta +\xi +\sigma a_h+{\gamma }_1+\mu },$

and offspring probability generating function (pgf) for resistant strain $z$, given $y(0)=0$ and $z(0)=1$ is

$f_2(u_1,u_2)=\frac{({\gamma }_2+\mu +d)+(\beta { }^{\mathrm{\prime }}+\xi )u_2^2}{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}.$

where $u_j,j=1,2$ denotes the individual in the compartment of sensitive strain and resistant strain, respectively. The power to $u_j$ denotes the number of offspring generated from one respective individual. Detailed calculations of the probability-generating functions are given in Appendix A2.

Remark 1. (1,1) is always a fixed point of $f_1$ and $f_2$, i.e.,

$f_1(1,1)=1$ and $f_2(1,1)=1$. Also, $f_i$‘s are not linear, i.e.,

$(1)f_i(u_1,u_2)\ne {\sum }_{j=1}^2{a}_j{u}_j;i=1,2,\dots .j=1,2,\dots$

$f_1(u_1,u_2)\ne a_1{u}_1+a_2{u}_2$

$f_2(u_1,u_2)\ne b_1{u}_1+b_2{u}_2$

$(2)f_i(0)\ne 0$

$f_1(0,0)=\frac{{\gamma }_1+\mu }{\beta +\xi +\sigma a_h+{\gamma }_1+\mu }\ne 0$

$f_2(0,0)=\frac{{\gamma }_2+\mu +d}{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}\ne 0$

Using the functions $f_1(u_1,u_2),f_2(u_1,u_2)$, we compute $n\times n$ non-negative matrix $\mathcal{M}=[m_{ij}{]}_{n\times n}$ called the Expectation Matrix. The entries of matrix $\mathcal{M}$ are obtained from partial differential equations given by:

$[m_{ij}{]}_{n\times n}=\frac{\mathrm{\partial }{f}_i}{\mathrm{\partial }{u}_j}{|}_{u=1}<\mathrm{\infty }$

This gives,

$\mathcal{M}=\left[\begin{array}{cc}\frac{2\beta }{\beta +\xi +\sigma a_h+\mu +{\gamma }_1}& 0\\ \frac{\xi +\sigma a_h}{\beta +\xi +\sigma a_h+\mu +{\gamma }_1}& \frac{2(\beta { }^{\mathrm{\prime }}+\xi )}{\beta { }^{\mathrm{\prime }}+\xi +\mu +{\gamma }_2+d}\end{array}\right]$

The expectation matrix $\mathcal{M}$ is irreducible [40,42]. Since off-spring pgfs are non-linear, and $\mathcal{M}$ is irreducible, there are at most two fixed points in [0,1] $\times$ [0,1], which depends on the spectral radius $\rho (\mathcal{M})$ of matrix $\mathcal{M}$. The branching process is called sub-critical, critical or super-critical if $\rho (\mathcal{M})<1$, $\rho (\mathcal{M})=1$ or $\rho (\mathcal{M})>1$, respectively. The existence of a non-trivial extinction probability in a parameter regime is determined by $\rho (\mathcal{M})>1$. We use the Jury condition to compute the parameter regime where $\rho (\mathcal{M})>1$ (see Appendix A3)

Now, suppose $\rho (\mathcal{M})>1$, and $(q_1,q_2)\in (0,1{)}^2$, we have a unique second fixed point. Then $f_2(q_1,q_2)=q_2$ gives

$\frac{({\gamma }_2+\mu +d)+(\beta { }^{\mathrm{\prime }}+\xi )q_2^2}{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}=q_2$

which produces quadratic equation in $q_2$

(9) $(\beta { }^{\mathrm{\prime }}+\xi )q_2^2-q_2(\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d)+({\gamma }_2+\mu +d)=0$(9)

Similarly $f_1(q_1,q_2)=q_1$ gives

$\frac{({\gamma }_1+\mu )+(\xi +\sigma a_h)q_2+\beta q_1^2}{\beta +\xi +\sigma a_h+{\gamma }_1+\mu +d}=q_1$

(10) $\beta q_1^2-(\beta +\xi +\sigma a_h+{\gamma }_1+\mu +d)q_1+({\gamma }_1+\mu )+(\xi +\sigma a_h)q_2=0$(10)

As colonization by drug-sensitive strain in the community inhibits invasion by drug-resistant [15], the probability of extinction of resistant strain with initial condition is given by

(11) $p=(1-q_1^{y(0)})q_2^{z(0)},0<q_1,q_2<1$(11)

where $q_1,q_2$ are obtained from equations (9)$(1)f_i(u_1,u_2)\ne {\sum }_{j=1}^2{a}_j{u}_j;i=1,2,\dots .j=1,2,\dots$ and (10$(2)f_i(0)\ne 0$).

5. Numerical computation and results

We implement the entire model given by equations (7)(7) $\frac{dh}{dt}=[r_h{h}^{\frac{\varsigma -1}{\varsigma }}+r_l{]}^{\frac{\varsigma }{\varsigma -1}}-{\delta }_{1}h$(7) as a stochastic system using a tau-leap algorithm, which is an event-driven approximated method of the Gillespie algorithm used for the simulation of stochastic systems [70]. In the process of developing antibiotic drug resistance, we consider a small population, where event-driven stochasticity is more relevant [71; 72]. The unit for the time step of the simulation is $day$. The algorithm updates the propensity functions at each interval length $\tau =day$ and then the reaction happens. Here, all demographic processes (i.e. birth and death) and disease processes (i.e. transmission, recovery, mutation and plasmid transfer), drug use, and growth of income in the deterministic model equations (7)(7) $\frac{dh}{dt}=[r_h{h}^{\frac{\varsigma -1}{\varsigma }}+r_l{]}^{\frac{\varsigma }{\varsigma -1}}-{\delta }_{1}h$(7) are assumed to be event-driven. Table 1 lists all the state transitions between model variables. The summary of tau-leap algorithm is as follows: let $g(t)$ be a disease compartment. Let $E_k$ be the list of events (such as the transmission of infection, recovery, mutation, plasmid transfer, etc.) that occurs at the rate $r_k(g(t))$ (most of our reactions are density dependent) and with a state change matrix $M_j$. First, we initialize the model with the initial condition $g(t_0)$ and then calculate the event rates $r_{k}g(t_0)$. We choose a fixed time step, i.e. $\tau =1$. For each event $E_k$, we then sample ${\xi }_k$ from a Poisson distribution with mean $\tau r_k$. This ${\xi }_k$ is the number of times each event occurs during the time interval $[t,t+\tau )$. Finally, we update the state by $g(t+\tau )=g(t)+{\sum }_k{\xi }_k{M}_j$. The initial populations for the stochastic run are considered as follows: $S(0)=20$, $y(0)=6$, $z(0)=2$, $h(0)=0.5$, with a total population of 200.

5.1. Mechanism of emergence and probability of extinction

We explore the extinction dynamics under two mechanisms: mutation and plasmid transfer. The mutation is omnipresent in the mechanism of evolution, but the ability of horizontal gene transfer differs tremendously between pathogens [78].

Figure 3(a) represents the time series of sensitive and resistant strains obtained from 10,000 sample runs of the model when both plasmid transfer and mutation drive the resistance separately, and simultaneously. The average population of resistant strain $z(t)$ grows and becomes extinct much faster in the case of plasmid transfer, which is not the same as the case of mutation. In the case of mutational resistance, both the growth and decay of resistance are relatively slow. The emergence of resistance through plasmid transfer follows a mass-action law, which ensures that the multiplication of resistance depends on the current densities of the resistant population $z(t)$ itself. So, the process develops faster. In contrast, mutation depends only on the density of sensitive strains $y(t)$ in the population. The continuous mutation from $y(t)$, though at a very small rate, ensures the longer-time persistence of resistant mutants in the community. We estimate the probability of extinction of resistant strains from these time-series simulations. For example, the mean time of extinction for Plasmid transfer is approximately 598 days, while for the case of mutation and both (Plasmid transfer & mutation) show approx. 700 days. Figure 3(b) also exhibits a box plot of this probability under different mechanisms over a longer period, and shows that it is much higher in the case of plasmid transfer. We observe a similar qualitative scenario when we consider the linear functional relation between economic status and antibiotic usage (Figure E3 in Appendix). However, the extinction of resistant strains provides a selective advantage for sensitive strain $y(t)$, but it also decays in the presence of antibiotic use. As the rate of overall plasmid transfer is relatively much higher than mutation, sensitive strain $y(t)$ also decays faster compared to the case of mutational resistance.

Figure 3. (a) Time series of sensitive strain $y(t)$ (Green) and resistant strain $z(t)$ (Red) and income $h(t)$ (Grey) under different mechanistic process of emergence: mutation, plasmid transfer, and both mutation and plasmid transfer. (b) Box plot of probability of extinction of resistant strain $z(t)$ under different process. We have simulated the model 5000 times for 700 days, and calculated the number of extinction at the ${700}^{th}$ day for three different processes with different initial conditions. In numerical case, we consider the extinction probability as proportion of runs which attains the value ‘0’ by 700 days.

5.2. Probability of disease extinction

Using parameter values given in Table 2, we numerically compute the probability of extinction of resistance strain under different initial values and compare the results obtained analytically from equation (11)(11) $p=(1-q_1^{y(0)})q_2^{z(0)},0<q_1,q_2<1$(11) . We consider 10,000 stochastic runs of the model under the same parameter regime with the same initial conditions and obtain that the proportion of resistant strain $z(t)$ hits zero. Table 3 presents the $p_0$, the probability of extinction under different initial values both analytically and numerically. As seen from the Table 3, both numerical approximations and analytical estimates of $p_o$ match with each other.

Table 2. Baseline parameter values for simulation.

Symbol	Description	Values (day${ }^{-1}$)	Refs
$\hat{\alpha }$	maximum antibiotic utilization rate	0.05	[73]
$\tilde{\alpha }$	minimum antibiotic utilization rate	0.001	[73]
$r_h^c$	intrinsic growth rate from capital	0.035	[61]
$r_l^c$	intrinsic growth rate from labour/technology	0.8	[61]
$h_o$	sensitivity of income to education and awareness	10–25	[59]
$p$	proportion of individuals who got recovered and take antibiotics	0.1	[74]
${\delta }_1$	rate of capital depreciation	0.0033	[61]
$1/{\delta }_2$	the average duration of antibiotic consumption	10 days	[75]
${\mathcal{R}}_{oy}$	reproduction number for sensitive strains	9.58	[75]
${\mathcal{R}}_{oz}$	reproduction number for drug-resistant strain	7.18	[75]
$1/{\gamma }_1$	the average duration of colonization of sensitive strain	60 days	[75]
$1/{\gamma }_2$	the average duration of colonization of resistant strain	90 days	[75]
$\mu$	natural mortality rate	0.00003629
$d$	death rate due to severity of resistant strain	0.0075	[75]
$\xi$	plasmid transfer rate	0.024	[76,77]
$\sigma$	mutation rate	${10}^{(-6)}$	[74,75]

Although we have not experimented with this extensively, the general pattern shows that the probability of extinction of resistance is higher when the initial number of the population with sensitive strain is high. This is also reflected in equation (11)(11) $p=(1-q_1^{y(0)})q_2^{z(0)},0<q_1,q_2<1$(11) . A probable ecological explanation for this is that the resistant strain cannot invade the population when the density of the sensitive strain is high. However, the probability of extinction is significantly low when the population with a resistant strain is initially greater than the population with a sensitive strain. The resistant pathogen circulates in the community and eventually establishes itself in the population by out-competing the sensitive strain.

5.3. Sensitivity analysis

A public health objective to eradicate or control the pathogen outbreak is to determine all pivotal characteristics of an epidemic that drives the infection. We compute partial rank correlation coefficients (PRCC) and perform a sensitivity analysis to determine the impact of uncertainty in the outcomes due to variations in certain parameters of the model [79–82]. PRCC measures the strength of the relationship between the model outcomes and the parameters, stating the degree of the effect that each parameter has on the outcome. We assume that all the model parameters are uniformly distributed. We have carried out 1000 simulations of the model using the ranges and baseline values given in Table 2.

The sensitivity analysis indicates the parameters that have the most influence on drug resistance transmission dynamics. Identification of given important parameters is necessary to form effective control measures for combating the spread of resistance. The parameter $h_o$ and ${\delta }_2$ are the two major parameters that influence the dynamics as seen by Figure 4.

Figure 4. Tornado plot (obtained using PRCC) exhibits influence of crucial parameters on the dynamics of variables $y(t)$, $z(t)$ and $h(t)$. The baseline parameters values are given in Table 2 and $h_o=1.5$, $\hat{a}=1.2$, $\tilde{a}=0.009$, $d=0$. For PRCC, ODE45 is used to analyse the solution for 0–700 days and the end point is taken to perform the sensitivity analysis.

The p-value of these two parameters helps us to understand the significance of the variables (Table 4).

Table 3. Probability of disease extinction $p$ from equation (11)(11) $p=(1-q_1^{y(0)})q_2^{z(0)},0<q_1,q_2<1$(11) and numerical approximation based on 10,000 sample time series of the model with parameter values $h_o=10$, $\sigma =0.00003$, $d=0.0098$.

$y(0)$	$z(0)$	Analytical ($p$)	Numerical
2	5	0.3138	0.3143
3	4	0.4062	0.4296
4	3	0.5113	0.5552
5	2	0.6400	0.6525
6	2	0.6402	0.6822

Table 4. $p$-value of determines the significance of parameters $h_o$ and ${\delta }_2$.

Parameter	$y(t)$	$z(t)$	$h(t)$
$h_o$	0.0026796	−0.022192	0.036996
${\delta }_2$	−0.0020894	0.033985	0.0125

Excessive use of antibiotics drives the growth of resistance [65]. The parameter ${\delta }_2$ reflects the average antibiotic use in the population, as $\frac{1}{{\delta }_2}$ is the mean duration of antibiotic use. As Figure 4(b) shows $z(t)$ is inversely proportional to ${\delta }_2$. Long-term use of antibiotics may fuel resistance in the community. The longer duration of antibiotic use gives rise to more $z(t)$ which declines the $y(t)$ and also income $h(t)$. ${\delta }_2$ is an important parameter in the dynamics that can be related to the public health control strategy. Restricting the duration of antibiotic use can slowdown the process of emergence of resistance by decreasing the selective pressure. In a study, it was observed that over-prescription of antibiotics was prevalent in patients for pulmonary infiltrates in ICU [83]. Drugs once started were continued for 4–20 days, disregarding the likelihood of pneumonia or even the absence of an extrapulmonary infection. A clinical control study over a small population in a hospital indicates that some patients diagnosed with recurrent pulmonary infections, who received antibiotic therapy only for 8 days show significantly less frequent emergence of multi-resistant pathogens as compared to patients who received treatment for the usual 20 days [83]. The parameter $h_o$ determines the impact of education and awareness of the community on antibiotic consumption. According to $a_h$ (Figure 2), the lower the value of $h_o$, the higher the volume of antibiotic use (see equation 4(4) $a_h(h)=(\hat{a}-\tilde{a})e^{-h{h}_o}+\tilde{a}$(4) ). Low education and awareness accelerate the imprudent practice of self-medication, which is more prevalent in developing countries. The likelihood of misuse is immense in low- and middle-income countries due to the high cost of treatments and also easy access to over-the-counter (OTC) drugs. Hence, $z(t)$ decreases with increase in $h_o$. A higher value of $h_o$ determines a low level of antibiotic use and low level of resistance, so $y(t)$ and $h(t)$ increase with the increase in value of $h_o$.

In particular, the results of this sensitivity analysis may provide insights into the public health authority's ability to design strategies that reduce the burden of drug resistance in the community. Public health programmes can use social media to spread awareness and educate people about the problem of drug resistance, put limitations on selling drugs over-the-counter and can provide incentives for proper treatment. Education and awareness can curtail significantly self-medication and over-the-counter (OTC) antibiotic consumption.

In Figure 4, ${\mathcal{R}}_{oy}$ is positively correlated with $y(t)$ and negatively correlated with $z(t)$ as more the transmission of sensitive strain, more is the increase in $y(t)$ population and decrease in the $z(t)$ population. ${\mathcal{R}}_{oz}$ is positively related with $z(t)$ as rising in the transmission of resistance strain will increase the population $z(t)$, which in turn, shrink the population of $y(t)$ and scale down the income $h(t)$ due to socio-economic impact of disease and infections. So, ${\mathcal{R}}_{oz}$ is negatively related to $y(t)$.

Resistance to antibiotics can occur either by mutations or by plasmid transfers, which are the two most important factors assumed in the current model of antibiotic drug resistance. With the increase in $\xi$ and $\sigma$, $z(t)$ increases (Figure 4(b)) and antibiotic consumption also increases while income decreases (Figure 4(c)). Antibiotics eliminate drug-sensitive strains, leaving behind resistant strains to proliferate due to natural selection [84]. So, $\xi$ and $\sigma$ have a negative impact on the population of $y(t)$ (Figure 4(a)).

5.4. Impact of ${\mathcal{R}}_{oy}$ and ${\mathcal{R}}_{oz}$ on probability of extinction

Basic reproduction numbers are very important parameters that impact the probability of extinction. The basic reproduction number of sensitive strain is given by:

${\mathcal{R}}_{oy}=\frac{\beta }{\mu +{\gamma }_1}$

The transmission potential of resistant strain is given by the basic reproduction number

${\mathcal{R}}_{oz}=\frac{\beta { }^{\mathrm{\prime }}}{\mu +p{\gamma }_2+d}$

We explore the possibility of changing the pattern in the probability of extinction depending on the basic reproduction ratio of two pathogens. Extensive and aberrant use of antibiotics exerts selective pressure, which provides sufficient conditions for a resistant strain to proliferate [85]. We shall assume that ${\mathcal{R}}_{oy}$ and ${\mathcal{R}}_{oz}$ are greater than one, but strains with a more basic reproduction number will dominate over the other strains.

From Figure 5, when ${\mathcal{R}}_{oy}>{\mathcal{R}}_{oz}$ sensitive strain will persist in the population. As a result, the population infected with resistance strain will be less, hence there will be a high probability of extinction of $z(t)$ (Figure 5). When ${\mathcal{R}}_{oy}<{\mathcal{R}}_{oz}$, the resistant strain easily dominates the sensitive strain, as the fitness cost for a resistant strain is less. The higher the value of ${\mathcal{R}}_{oz}$, the easier the elimination of sensitive strain. The community will be invaded by the resistant strain. With values ${\mathcal{R}}_{oz}\in (3-4)$, both sensitive and resistant strains coexist in the population – this region determines the medium probability of extinction. Thus, for high extinction of $z(t)$ basic reproduction number of the sensitive strain should be greater than the reproduction number of the resistant strain to existing in the population in the presence of the resistance strain. Together, it is observed that ${\mathcal{R}}_{oz}$ is more influential in changing the extinction dynamics compared to ${\mathcal{R}}_{oy}$. However, this pattern also depends on other socio-economic parameters such as poverty, income, self-medication, poor community hygiene, and so on.

Figure 5. Impact of ${\mathcal{R}}_{oy}$ and ${\mathcal{R}}_{oz}$ on extinction of resistance strain given by the contour plot. It depicts that for lower value of ${\mathcal{R}}_{oz}$ only sensitive strain exist so probability of extinction is high. When fitness of resistance strain increases ${\mathcal{R}}_{oz}=6-12.5$, strains starts multiplying and thus invade the sensitive strain.

5.5. Impact of the volume of antibiotic consumption, awareness, and education level $h_o$ on the probability of extinction

Although sensitivity analysis indicates a high correlation of drug resistance with the volume of drug use, in this section, we quantify the nonlinear dependence of these variates by portraying the probability of extinction of resistant strain in the two-parameter regime, the awareness level $h_o$ and the variability of drug use $(\hat{a}-\tilde{a})$ by the community.

Figure 6 shows the contour plot of probabilities of extinction for $h_o$ and difference $(\hat{a}-\tilde{a})$ of antibiotic consumptions. The image follows a pattern from the lower to the higher value of the probability of extinction of $z(t)$. It exhibits a high probability of extinction when the difference $(\hat{a}-\tilde{a})$ is low, but there is a little impact of $h_o$. In contrast, there is a high non-linear dependence of $h_o$, when the difference is relatively higher. The probability of extinction increases from low to moderate as the $h_o$ increases when the difference in antibiotic use is high. With higher $h_o$ the total volume of antibiotic use decreases, and hence the fitness level of $z(t)$ shrinks ensuring the co-existence of strains, and so, the probability of extinction is medium.

Figure 6. The contour plot of resistance dynamics depending the differences in maximum and minimum drug usages $\hat{a}(=a_{max})-\tilde{a}(=a_{min})$ and $h_o$. The low education in developing nations drives higher antibiotic use. In developed nations difference in drug usage is low, and hence $h_o$ has negligible impact. The increased difference $\hat{a}-\tilde{a}$ implies a low-middle income countries where education and knowledge about antibiotic use plays an important role. For higher values in difference and $h_o=10-15$, density of population infected with resistance strain is more, which makes the extinction difficult. Baseline parameters are given in Table 2.

We may characterize this non-linear pattern of the probability of extinction in terms of the development profile of nations. As indicated by the red star and white circle, we can identify that developing countries are apparently situated in the left-uppermost corner of the diagram, whereas developed populations are placed in the right-lowermost corner. Extinction is more unlikely in developing nations. Partly, because of the low economy and inappropriate use of antibiotics increase the burden of disease in such communities, which also inhibits socio-economic growth and further promotes self-medication in the population. In fact, Klein et al. (2018) indicates in a study that the daily dosage per person has increased by 65% in India and 70% in China since 2000. If other countries follow this path, the researchers pointed out that global antibiotic use could reach 125 billion doses per year, up from 42 billion in 2015. In a recent study, antibiotic use increased 65%, and antibiotic consumption rate increased 39% [17]. Thus, education and awareness about self-medication play an important factor in the establishment and removal of resistant strains from the community, and the public health community should focus on campaigning more awareness programmes in developing populations.

6. Discussion

The biological process of drug resistance in pathogens can be divided into two phases: i) the individual-level process is the emergence or appearance of resistant mutants in healthy individuals, and ii) the population-level phenomenon – the establishment of the pathogen [78]. Thus, public health intervention should be a combination of both individual and population level.

The proliferation of resistance and its abundance in society is closely dependent on socioeconomic status, which is a composite measure that typically incorporates economic status (measured by income), social status (measured by awareness and education), and work status (measured by occupation) [86]. In the current analysis, we hypothesize an exponential function to represent the reciprocal relation between socioeconomic status and self-medication in the population, and examine the impact of such factors on the emergence and extinction of resistance by developing a non-linear CTMC model integrating socioeconomic growth and antibiotic misuse. A multi-type branching process is considered to approximate the nonlinear CTMC and estimate the probability of disease outbreak or extinction (Table 3).

Both analytical derivation and numerical computation of the probabilities of extinction of resistance strains mirror each other for different initial sizes of infectives in the population. The results clearly indicate that the probability of eliminating resistance is high if initially fewer people are incurred with resistance as opposed to the situation where , initially, a large number of people are already infected with resistance strain. It will be difficult to control the situation, once resistance gets established in the population. The analysis also reveals that interventions to control the spread of resistance should be initiated at the beginning of the outbreak. The sensitivity analysis strongly depicts that public health authorities should focus on improvements in the cost of treatment and health facilities, putting control on overuse and misuse of antibiotics including OTC drug sales and increasing awareness programmes of antibiotic resistance among the common public. Knowledge and delving into the biology of antibiotic drug resistance have a significant impact on the implementation and development of policies regarding infection control [87–89].

Estimating the epidemic threshold for the occurrence of a disease outbreak or extinction is an important issue in public health. This is especially true for drug resistance model integrating socioeconomic factors. Every country has its own economic constraints and epidemiological parameters behind the emergence and extinction of drug resistance in its population. Also, different pathogens may have different mechanisms of emergence and selection. Despite these, our stochastic model provides a framework to analyse the likelihood of extinction or establishment of a resistant pathogen depending on income, amount of drug use, epidemiological factors, and others. The PRCC analysis determines parameter dependencies on the probability of extinction of the resistant pathogens. In reality, the association between irrational antibiotic consumption and the socioeconomic growth of the population is more complex and nonlinear, and there is rarely any empirical contribution in support of this. So, we considered a simpler approach to illustrate the complex outcomes as emergent properties of the integrated model of emergence and establishment of drug resistance and socioeconomic factors of the population.

There are scopes to improve the current model and analysis. In reality, resistance may emerge and spread through other means like hospital infection or due to unhygienic conditions. Here, we have only considered antibiotic misuse. A detailed model may be developed to address questions on the emergence and extinction of drug resistance, possibly integrating individual and population-level dynamics of the emergence and establishment of the resistant pathogens. The reciprocal relation between income and antibiotic misuse may not always be true for any community. There are many socio-economic and policy-oriented factors that define the complex relation. One can use the same model framework and investigate the model outcome under different nonlinear relations of income and antibiotic usage through self-medications. Also, the incorporation of health risk behaviours, migration, and case importation in the model framework might increase the accuracy of predictions, which would provide insights into public health for combating resistance globally.

Acknowledgments

The authors thank the Science and Engineering Research Board (SERB), DST, India for financial support (ECR/2016/000591) to conduct this research. We thank Associated Editor and reviewers of the journal for useful comments to represent the materials better in the draft. We also thank the DST-FIST lab in SNIoE for facilitating the computational facilities required for this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is funded by the Science and Engineering Research Board (SERB), DST, India (ECR/2016/000591).

Appendix A

A1. Next-generation matrix

The basic reproductive number, ${\mathcal{R}}_o$ is the most important metric in disease modelling, and one of the most valuable ideas that mathematical thinking has brought to epidemic theory [67,90]. The magnitude of ${\mathcal{R}}_o$ is used to measure the risk of an epidemic in emerging infectious diseases. It is defined as the expected number of secondary cases produced by a single infection in a totally susceptible population. The next-generation approach is a more general method of computing ${\mathcal{R}}_o$ for SIR-type extended models [64–66]. In next-generation method, ${\mathcal{R}}_o$ is defined as the spectral radius of the ‘next-generation operator’ [67].

Consider there are $n$ compartments out of which $m$ are infected. Define the vector $\bar{x}=x_i,i=1,..,n$, where $x_i$ denotes the number or proportion of individuals in the $ith$ compartment. Let $F_i(\bar{x})$ be the rate of appearance of new infections in compartment $i$ and let $V_i(\bar{x})=V_i^{-}(\bar{x})-V_i^{+}(\bar{x})$, where $V_i^{+}$ is the rate of transfer of individuals into compartment $i$ by all other means and $V_i^{-}$ is the rate of transfer of individuals out of $i^{th}$ compartment. The difference $F_i(\bar{x})-V_i^{-}(\bar{x})$, gives the rate of change of $x_i$. $F_i$ include infections that are newly arising, but does not include terms which describe transfer of infectious individuals from one infected compartment to other. $V_i$ denotes transfers of infections from one compartment to another. We can form the next generation matrix $F{V}^{-1}$ from matrices of partial derivatives of $F_i$ and $V_i$.

$F=\frac{\mathrm{\partial }{F}_i(x_o)}{\mathrm{\partial }{x}_j}$

and

$V=\frac{\mathrm{\partial }{F}_i(x_o)}{\mathrm{\partial }{x}_j}$

where $x_o$ is the disease-free equilibrium. The entries of $F{V}^{-1}$ give the rate at which infected individuals in $x_j$ produce new infections in $x_i$, times the average length of time an individual spends in a single visit to compartment $j$. ${\mathcal{R}}_o$ is given by the spectral radius (dominant eigenvalue) of the matrix $F{V}^{-1}$ [91–94].

A2. Computation of pgfs for the drug-resistance model

Let time be continuous, $t\in [0;\mathrm{\infty })$ and let $S(t)$, $y(t)$, $z(t)$, $h(t)$, represent discrete random variables for the number of susceptibles with antibiotics, the number of infected individuals with sensitive strain and resistant strain, respectively. We consider a small population, such that all variables are bounded.

Since $y(t)$ and $z(t)$ are the only variables that are directly responsible for the emergence of resistance in the community, we apply the multi-type branching process and compute pgfs for these two variables.

The offspring pgf for $y(t)$ given $y(0)=1$ and $z(t)=0$ is given by:

$f_1(u_1,u_2)=Prob\{y_{11}=k_1,y_{12}=k_2\}$

$=\sum _{k_2=0}^{\mathrm{\infty }}\sum _{k_1=0}^{\mathrm{\infty }}Prob\{k_1,k_2\}{u}_1^{k_1}{u}_2^{k_2}$

$=Prob\{k_1=0,k_2=0\}{u}_1^0{u}_2^0+Prob\{k_1=0,k_2=1\}{u}_1^0{u}_2^1+$

$Prob\{k_1=1,k_2=0\}{u}_1^1{u}_2^0+Prob\{k_1=1,k_2=1\}{u}_1^1{u}_2^1+$

$Prob\{k_1=2,k_2=0\}{u}_1^2{u}_2^0$

$f_1(u_1,u_2)=\frac{({\gamma }_1+\mu )+(\xi +\sigma a_h)u_2+\beta u_1^2}{\beta +\xi +\sigma a_h+{\gamma }_1+\mu }$

To describe the probabilities in the calculation above, one individual with sensitive strain dies or is removed with probability $\frac{({\gamma }_1+\mu )}{\beta +\xi +\sigma a_h+{\gamma }_1+\mu }$, this means an individual with sensitive strain dies or is removed before infecting with resistance strain. The $y(t)$ giving birth to $y(t)$ with probability $\frac{\beta }{\beta +\xi +\sigma a_h+{\gamma }_1+\mu }$. This means that an individual with a sensitive strain infects another individual with a sensitive strain which results in $(X_{11}=2,X_{12}=0)$. The term $u_1^2$ means an individual with a sensitive strain produces two individuals with sensitive strains. A single $y(t)$ gives birth to a single individual with resistance with probability $\frac{(\xi +\sigma a_h)}{\beta +\xi +\sigma a_h+{\gamma }_1+\mu }$.

Similarly, the pgf for $z(t)$, given $y(0)=0$ and $z(t)=1$

$f_2(u_1,u_2)=Prob\{y_{21}=k_1,y_{22}=k_2\}$

$=\sum _{k_2=0}^{\mathrm{\infty }}\sum _{k_1=0}^{\mathrm{\infty }}Prob\{k_1,k_2\}{u}_1^{k_1}{u}_2^{k_2}$

$=Prob\{k_1=0,k_2=0\}{u}_1^0{u}_2^0+Prob\{k_1=0,k_2=1\}{u}_1^0{u}_2^1+$

$Prob\{k_1=1,k_2=0\}{u}_1^1{u}_2^0+Prob\{k_1=1,k_2=1\}{u}_1^1{u}_2^1+$

$Prob\{k_1=0,k_2=2\}{u}_1^2{u}_2^0$

$f_2(u_1,u_2)=\frac{({\gamma }_2+\mu +d)+(\beta { }^{\mathrm{\prime }}+\xi )u_2^2}{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}.$

For $f_2$, one individual with resistant strain dies or is removed with probability $\frac{({\gamma }_2+\mu +d)}{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}$. An individual with a resistant strain infects another individual with the resistant strain which results in $(X_{21}=0,X_{12}=0$ with probability $\frac{\beta { }^{\mathrm{\prime }}+\xi }{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}$. The term $u_2^2$ means an individual with resistant strain produces two individuals.

A3. Jury condition and computation of extinction probability

Proposition.

A non-negative square matrix A of order $n$ is irreducible $\iff$ $(I+A{)}^{n-1}>0$, where I is n $\times$ n identity matrix and $\prime$ stands element wise larger than.

Since the off-spring pgf is non-singular or not simple and the expectation matrix $\mathcal{M}$ is irreducible, then there are at most two fixed points in [0,1] $\times$ [0,1]. The probability of disease extinction or persistence depends on the size of the spectral radius of matrix $\mathcal{M}$, $\rho (\mathcal{M})$. The branching process is called subcritical, critical or supercritical if $\rho (\mathcal{M})<1$, $\rho (\mathcal{M})=1$ or $\rho (\mathcal{M}>)1$ respectively. We are interested in a parameter regime of our model, where it is supercritical and we have nontrivial probabilities of disease extinction. Below, we use the Jury condition to determine the supercritical parameter regime.

(1) If $\rho (M)<1$ or $=1$ then (1,1) is only fixed point of the system and the probability of disease extinction is 1.

Jury conditions: The necessary and sufficient condition for a real 2 $\times$ 2 matrix A to have eigenvalue with modulus less than 1 is

$Tr(A)<1+Det(A)<2$

The 2$\times$ 2 matrix $\mathcal{M}$ satisfies $\rho (\mathcal{M})<1\iff$ jury conditions hold.

Now,

$Tr(\mathcal{M})=\frac{2\beta }{\beta +\xi +\sigma a_h+{\gamma }_1+\mu }+\frac{2(\beta { }^{\mathrm{\prime }}+\xi )}{\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d}$

$Det(\mathcal{M})=\frac{4\beta (\beta { }^{\mathrm{\prime }}+\xi )}{(\beta +\xi +\sigma a_h+{\gamma }_1+\mu )(\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d)}$

(a) The inequality $1+Det(\mathcal{M})<2Det\mathcal{M}<1$, which gives

(12) $4\beta (\beta { }^{\mathrm{\prime }}+\xi )<AB$(12)

(b) The inequality $Tr(\mathcal{M})<1+Det(\mathcal{M})$ implies

$\frac{2\beta }{A}+\frac{2(\beta { }^{\mathrm{\prime }}+\xi )}{B}<1+\frac{4\beta (\beta { }^{\mathrm{\prime }}+\xi )}{AB}$

$\Rightarrow 2\beta B+2(\beta { }^{\mathrm{\prime }}+\xi )A<AB+4\beta (\beta { }^{\mathrm{\prime }}+\xi ),$

which finally shows,

(13) $\frac{\beta }{A}+\frac{(\beta { }^{\mathrm{\prime }}+\xi )}{B}<1.$(13)

So, assuming our parameters are in the complement of this regime given by equations (12)(12) $4\beta (\beta { }^{\mathrm{\prime }}+\xi )<AB$(12) and (13(13) $\frac{\beta }{A}+\frac{(\beta { }^{\mathrm{\prime }}+\xi )}{B}<1.$(13) ), we have the following:

(2) $\rho (\mathcal{M})>1$ and there exists a unique second fixed point $(q_1,q_2)\in (0,1{)}^2$, which determines the existence or persistence of disease.

Considering $f_2(q_1,q_2)=q_2$ gives

$({\gamma }_2+\mu +d)+(\beta { }^{\mathrm{\prime }}+\xi )q_2^2=q_2(\beta { }^{\mathrm{\prime }}+\xi +{\gamma }_2+\mu +d)$

Similarly, $f_1(q_1,q_2)=q_1$ gives

$\beta q_1^2-(\beta +\xi +\sigma a_h+{\gamma }_1+\mu +d)q_1+({\gamma }_1+\mu )+(\xi +\sigma a_h)q_2=0$

Thus, the probability of extinction of the resistant strain with the initial condition is given by:

$p=(1-q_1^{y(0)})q_2^{z(0)}.$

Figure E1. Plots of the antibiotic use (left panel) and resistance prevalence (right panel) for different antibiotics. These data are for some of the representative developing nations which compares the antibiotic use and level of resistance corresponding to certain selected antibiotics. The time series shows that the increment in antibiotics use over the years increases continuously, which increases the resistance level in these nations over the years.

Figure E2. SMA prevalence in different countries: India (blue, 1–7 [22–27]), Indonesia (cyan, 8–10 [28–30]), Korea (green, 11 [31]), Nepal (grey, 12–14 [32–34]), Sri Lanka (black, 15 [35]), Thailand (maroon, 16 [36]), Bangladesh (purple, 17–18 [37,38]), Bhutan (yellow, 19 [39]). Data obtained from [21].

Figure E3. (a) Time series of sensitive strain $y(t)$ (Green) and resistant strain $z(t)$ (Red) and income $h(t)$ (Grey) under different mechanistic process of emergence: mutation, plasmid transfer, and both mutation and plasmid transfer. (b) Box plot of probability of extinction of resistant strain $z(t)$ under different process.