Mathematical modelling of antibiotic interaction on evolution of antibiotic resistance: an analytical approach

Model formulation

Portions of this text were previously published as part of a preprint (Nashebi, Sari & Kotil, 2022). In our study, we focus on modeling the scenario of multidrug treatment against the Staphylococcus aureus bacteria in an individual. The population sizes of wildtype and mutant bacteria at time t are denoted as S(t) and R(t), respectively. We assume that bacterial growth follows a logistic model with a carrying capacity K. The birth rate of wildtype bacteria is represented by β_s, and the birth rate of mutant bacteria is denoted as β_r. It is important to note that specific mutations that confer resistance to chemical control incur a fitness cost, which can result in reduced reproductive capacity or competitive ability (Alavez-Ramírez et al., 2007). We quantify this fitness cost as a reduction in the reproduction rate of the mutant strain, leading to β_r ≤ β_s. The natural death rates for wildtype and mutant bacteria are represented as μ_s and μ_r, respectively. Additionally, both bacterial types can also die due to the action of antibiotics. In this study, we consider two types of antibiotics: (1) antibiotic agent M, which kills both wildtype and mutant bacteria, (2) antibiotic agent N, which also kills both wildtype and mutant bacteria. The antibiotic M kills wildtype and mutant bacteria with rates ${\alpha }_{11}$ and ${\alpha }_{21}$, while the antibiotic N affects them with rates ${\alpha }_{12}$ and ${\alpha }_{22}$, respectively. Moreover, these antibiotics interact synergistically and antagonistically to kill wildtype and mutant bacteria. To account for the combined effect of antibiotics M and N in killing both wildtype and mutant bacteria, we adopt a density-dependent approach based on the work of Ibargüen-Mondragón et al. (2014). This allows us to explore the relationship between the pharmacodynamics of the antibiotics and the population dynamics of wildtype and mutant bacteria when exposed to these antibiotics. The model is described as follows:

(1a) $\overline{X_s}=\left(\overline{{\alpha }_{11}}{C}_1+\overline{{\alpha }_{12}}{C}_2+{\lambda }_1\overline{{\alpha }_{11}}\overline{{\alpha }_{12}}{C}_1{C}_2\right)$

(1b) $\overline{X_r}=\left(\overline{{\alpha }_{21}}{C}_1+\overline{{\alpha }_{22}}{C}_2+{\lambda }_2\overline{{\alpha }_{21}}\overline{{\alpha }_{22}}{C}_1{C}_2\right)$where

(2a) $\overline{{\alpha }_{11}}={\displaystyle \frac{E_{max}^{M,S}}{I{C}_{50}^{M,S}}}$

(2b) $\overline{{\alpha }_{12}}={\displaystyle \frac{E_{max}^{N,S}}{I{C}_{50}^{N,S}}}$and

(3a) $\overline{{\alpha }_{21}}={\displaystyle \frac{E_{max}^{M,R}}{I{C}_{50}^{M,R}}}$

(3b) $\overline{{\alpha }_{22}}={\displaystyle \frac{E_{max}^{N,R}}{I{C}_{50}^{N,R}}}$where, E^M,S_max and E^N,S_max represent the maximal killing rates of antibiotics M and N of wildtype bacteria, E^M,R_max and E^N,R_max represent the maximal killing rates of antibiotics M and N of mutant bacteria. IC^M,S₅₀ and IC^N,S₅₀ signify the half-maximal inhibitory concentration of the antibiotics $M$ and N for wildtype bacteria. IC^M,S₅₀ and IC^N,S₅₀ denote the half-maximal inhibitory concentration of the antibiotics M and N for mutant bacteria. Multiplication of Eqs. (2a), (2b) and (3a), (3b) give Emax model (Salahudeen & Nishtala, 2017; Holford, 2017) which quantify the relationship between concentration and effect of antibiotics.

The parameters λ₁ and λ₂ represent the interaction strengths for wildtype and mutant bacteria, respectively. These parameters have a range of −1.5 to 1.5 (−1.5 ≤ λ₁, λ₂ ≤ 1.5) (Torella et al., 2010). Negative values indicate antagonistic interactions between the antibiotics, while positive values indicate synergistic interactions.

Bacteria have the potential to acquire resistance to both antibiotic agents through mutation. The concentrations of antibiotics M and N are denoted as C₁(t) and C₂(t), respectively. We assume that these two antibiotics belong to the same class and have the same inhibitory effect. Furthermore, we assume that both antibiotics bind to the same target, allowing bacteria to develop resistance to both antibiotics through a single mutation event. The acquisition of resistance by mutant bacteria from wildtype bacteria is modeled by the terms q₁C₁(t)S(t) and q₂C₂(t)S(t), where q₁ and q₂ represent the mutation rates of wildtype bacteria when exposed to antibiotics, respectively.

To maintain a constant concentration of antibiotics M and N, they are supplied at a constant rate θ₁ and θ₂, respectively. Antibiotics are removed from the system at a constant per capita rate μ₁ and μ₂, respectively.

Under the assumptions revealed above, we obtain the following system of differential equations:

(4a) ${\displaystyle \frac{dS}{dt}={\beta }_{s}S\left(1-{\displaystyle \frac{S+R}{K}}\right)-\left(q_1{C}_1+q_2{C}_2\right)S-\left(\overline{X_s}+{\mu }_s\right)S}$

(4b) ${\displaystyle \frac{dR}{dt}={\beta }_{r}R\left(1-{\displaystyle \frac{S+R}{K}}\right)+\left(q_1{C}_1+q_2{C}_2\right)S-\left(\overline{X_r}+{\mu }_r\right)R}$

(4c) ${\displaystyle \frac{d{C}_1}{dt}={\theta }_1-{\mu }_1{C}_1}$

(4d) ${\displaystyle \frac{d{C}_2}{dt}={\theta }_2-{\mu }_2{C}_2}$with the consideration of following change of variable:

$s={\displaystyle \frac{S}{K},r=\frac{R}{K},c_1=\frac{C_1}{{\theta }_1/{\mu }_1},c_2={\displaystyle \frac{C_2}{{\theta }_2/{\mu }_2}}}$the non-dimensionalized system (4) can be rewritten as:

(5a) ${\displaystyle \frac{ds}{dt}={\beta }_{s}s\left(1-\left(s+r\right)\right)-(q_1{c}_1+q_2{c}_2)s-\left(X_s+{\mu }_s\right)s}$

(5b) ${\displaystyle \frac{dr}{dt}={\beta }_{r}r\left(1-\left(s+r\right)\right)+(q_1{c}_1+q_2{c}_2)s-\left(X_r+{\mu }_r\right)r}$

(5c) ${\displaystyle \frac{d{c}_1}{dt}={\mu }_1-{\mu }_1{c}_1}$

(5d) ${\displaystyle \frac{d{c}_2}{dt}={\mu }_2-{\mu }_2{c}_2}$where

(6a) $X_s=\left({\alpha }_{11}{c}_1+{\alpha }_{12}{c}_2+{\lambda }_1{\alpha }_{11}{\alpha }_{12}{c}_1{c}_2\right)$

(6b) $X_r=\left({\alpha }_{21}{c}_1+{\alpha }_{22}{c}_2+{\lambda }_2{\alpha }_{21}{\alpha }_{22}{c}_1{c}_2\right).$and

$${\alpha }_{1i}=\overline{{\alpha }_{1i}}\left({\theta }_i/{\mu }_i\right)$$

${\alpha }_{2i}=\overline{{\alpha }_{2i}}\left({\theta }_i/{\mu }_i\right),i=1,2$

The region of biological interest of system (5) is given by

(7) $\mathrm{\Omega }=\left\{\left(s,r,c_1,c_2\right)ϵR_{+}^4:0\le s,r,c_1,c_2\le 1,0\le s+r\le 1\right\}.$

The set Ω defined in (7) is positively invariant for the system (5) (Ibargüen-Mondragón et al., 2019). Consequently, the system (5) is well-posed because solutions with initial conditions in Ω remain there for all t ≥ 0.

Qualitative analysis of the model

In this part, we will analyze the solutions of system (5) which include infection-free, all-mutant, and coexistence of wildtype and mutant bacteria equilibrium-points. We will then investigate the stability conditions of these solutions based on the antibiotic interaction parameters (λ1 and λ2) for wildtype and mutant bacteria.

Equilibrium solutions

The model (5) always contains the infection-free equilibrium P₀ = (0,0,1,1) in Ω. This equilibrium point represent state where both wildtype and mutant bacteria are eliminated under combination therapy. If R_r > 1, P₁ = (0, (R_r − 1)/R_r,1,1) is an all-mutant equilibrium in Ω where

(8) $R_r={\displaystyle \frac{{\beta }_r}{\left({\alpha }_{21}{\alpha }_{22}{\lambda }_2+{\alpha }_{21}+{\alpha }_{22}\right)+{\mu }_r}.}$

This equilibrium point represent state where only mutant bacteria persist under combination therapy. When R_s > 1 and R_s > R_r in addition to P₀, and P₁ there exists a coexistence of wildtype and mutant bacteria equilibrium in Ω, P₂ ( $\overline{s}$, $\overline{r}$,1,1) where

(9) $R_s={\displaystyle \frac{{\beta }_s}{m+\left({\alpha }_{11}+{\alpha }_{12}+{\lambda }_1{\alpha }_{11}{\alpha }_{12}\right)+{\mu }_s},}$

(10) $\overline{r}={\displaystyle \frac{m\left({\displaystyle \frac{R_s-1}{R_s}}\right)}{{\beta }_r\left({\displaystyle \frac{1}{R_r}-{\displaystyle \frac{1}{R_s}}}\right)+m},}$and

(11) $\overline{s}={\displaystyle \frac{R_s-1}{R_s}-r.}$

This equilibrium point represent state where both wildtype and mutant bacteria persist under combination therapy. The derivation of equilibrium points has been given in Supplemental Information.

Based on the traditional definition of the basic reproduction number, this quantity

(12) $N_r={\displaystyle \frac{{\beta }_r}{{\mu }_r},}$

This parameter is construed as the product of the mutant bacteria’s reproduction rate (β_r) and their average lifespan (1/μ_r). It signifies the count of bacteria generated by a mutant bacterium throughout its typical lifetime. Likewise,

(13) $$N_s={\displaystyle \frac{{\beta }_s}{{\mu }_s}}$$

It is understood as the quantity of bacteria generated by a wildtype bacterium during its average lifespan. Conversely, Rs, as defined in Eq. (9), is redefined as

(14) $R_s={\displaystyle \frac{{\mu }_s}{m+\left({\alpha }_{11}+{\alpha }_{12}+{\lambda }_1{\alpha }_{11}{\alpha }_{12}\right)+{\mu }_s}{N}_s.}$

Since

(15) $${\displaystyle \frac{{\mu }_s}{m+\left({\alpha }_{11}+{\alpha }_{12}+{\lambda }_1{\alpha }_{11}{\alpha }_{12}\right)+{\mu }_s}}$$specify the proportion of wildtype bacteria that haven’t undergone spontaneous mutations and remain unaffected by antibiotics. Therefore, R_s represents the count of bacteria produced by this fraction of wildtype bacteria that hasn’t undergone spontaneous mutations and remains unaffected by the combination therapy. Similarly, R_r as defined in Eq. (8) is reformulated as

(16) $$R_r={\displaystyle \frac{{\mu }_r}{\left({\alpha }_{21}+{\alpha }_{22}+{\alpha }_{21}{\alpha }_{22}{\lambda }_2\right)+{\mu }_r}{N}_r}$$

This represents the count of bacteria produced by the fraction of mutant bacteria that evade the effects of combination therapy.

The results demonstrate the following: (a) when the average bacteria count produced by the fraction of mutant bacteria evading the antibiotic combination effect exceeds one (R_r > 1), the population of mutant bacteria will endure, (b) if the average bacteria count produced by the fraction of wildtype bacteria without mutations, evading the antibiotics’ combination effect, is greater than one (R_s > 1), and the average bacteria count produced by the fraction of mutant bacteria exceeds one, both susceptible and resistant bacteria will persist.

Stability of equilibria points

In this section, we determine the local asymptotic stability of the equilibrium solutions of the system (5). Linearization of the system (5) around point P is given by:

(17) $${\overrightarrow{x}}^{\mathrm{\prime }}=J\left(P\right)\overrightarrow{x}$$where

(18) $$\overrightarrow{x}={\left(s,r,c_1,c_2\right)}^T$$and the matrix J evaluated at P is:

(19) $$J\left(P\right)=\left[\begin{array}{cccc}{j}_{11}\left(P\right)& -{\beta }_{s}s& -\left(c_2{\lambda }_1{\alpha }_{11}{\alpha }_{12}+{\alpha }_{11}\right)s& -\left(c_1{\lambda }_1{\alpha }_{11}{\alpha }_{12}+{\alpha }_{12}\right)s\\ -{\beta }_{r}r+m& j_{22}\left(P\right)& -\left(c_2{\lambda }_2{\alpha }_{21}{\alpha }_{22}+{\alpha }_{21}\right)r& -\left(c_1{\lambda }_2{\alpha }_{21}{\alpha }_{22}+{\alpha }_{22}\right)r\\ 0& 0& -{\mu }_1& 0\\ 0& 0& 0& -{\mu }_2\end{array}\right]$$with

(20a) $j_{11}\left(p\right)={\beta }_s\left(1-\left(s+r\right)\right)-{\beta }_{s}s-m-\left(\left({\alpha }_{11}{c}_1+{\alpha }_{12}{c}_2+{\lambda }_1{\alpha }_{11}{\alpha }_{12}{c}_1{c}_2\right)+{\mu }_s\right)$

(20b) $j_{22}\left(p\right)={\beta }_r\left(1-\left(s+r\right)\right)-{\beta }_{r}r-\left(\left({\alpha }_{21}{c}_1+{\alpha }_{22}{c}_2+{\lambda }_2{\alpha }_{21}{\alpha }_{22}{c}_1{c}_2\right)+{\mu }_r\right).$

By evaluating the Eq. (19) Jacobian J in P₀, P₁ and P₂ (see Supplemental Information) we obtain that, firstly, if R_s < 1 and R_r < 1, then the infection-free equilibrium P₀ is locally and asymptotically stable in Ω. If R_s > 1 or R_r > 1, then P₀ is unstable. Since α₁₁, α₁₂, μ_s, and β_s are positive; there are three conditions for R_s < 1 if λ₁ > 0, λ₁ < 0, or λ₁ = 0. If λ₁ > 0, λ₁ < 0 the necessary condition for R_s < 1 is:

(21) $${\beta }_s-{\mu }_s-m<{\alpha }_{11}+{\alpha }_{12}+{\lambda }_1{\alpha }_{11}{\alpha }_{12}$$and if λ₁ = 0, the necessary condition is:

(22) ${\beta }_s-{\mu }_s-m<{\alpha }_{11}+{\alpha }_{12}.$

This implies that when antibiotics combination eliminates the wildtype bacteria, and prohibit the proliferation of mutants, in this case, both bacteria die out. Secondly, If R_r > R_s and R_r > 1, then the equilibrium P₁ is locally and asymptotically stable in Ω. If R_r < R_s or R_r <1, then P₁ is unstable. Since α₂₁, α₂₂, μ_r, and β_r are positive, there are three conditions for R_r > 1, if λ₂ > 0, λ₂ < 0, or λ₂ = 0. If λ₂ > 0, λ₂ < 0 the necessary condition for R_r < 1 is:

(23) $${\beta }_r-{\mu }_r>{\alpha }_{21}+{\alpha }_{22}+{\lambda }_2{\alpha }_{21}{\alpha }_{22}$$and if λ₂ = 0 the necessary condition is:

(24) ${\beta }_r-{\mu }_r>{\alpha }_{21}+{\alpha }_{22}.$

In this scenario, assuming mutants have an average reproduction rate greater than one and the reproductive capacity of wildtype bacteria is lower than that of mutants, only mutants survive while wildtype bacteria go extinct.

Finally, if R_s > 1 and R_s > R_r, the equilibrium point P₂ is within Ω and is both locally and asymptotically stable. Here, when wildtype bacteria have a reproduction rate greater than one and a higher reproductive capacity than mutants, both strains can coexist. Despite the lower reproductive capacity of mutants compared to wildtype bacteria, the occurrence of spontaneous mutations in wildtype strains enables their survival.

Numerical simulations

This section gives some numerical justification for the equilibrium points and their stability criterion. Since these equilibrium points demonstrate the free-infection, all-mutant, and coexistent of wildtype and mutant bacteria conditions for bacteria population, it is important to have some numerical justification for them. The parameters used in the simulations are constant and are given in Table 1. For the numerical simulation, we consider an individual with a disease caused by Staphylococcus aureus bacteria that develop resistance to antibiotics M and N through mutation. The antibiotic interaction parameter for wildtype and mutant bacteria are (λ₁) and (λ₂), respectively. As underlined in the work of Torella et al. (2010), λ equals $0$ for additive interaction, 1 for synergistic interaction, and −1 for antagonistic interaction. Our simulation follows three scenarios. In the first scenario, antibiotics interact additively against wildtype and the mutant bacteria (λ₁ = λ₂ = 0). In the second scenario, antibiotics interact synergistically with the wildtype bacteria but interact antagonistically with mutants (λ₁ = 1, λ₂ = −1). In the third scenario, antibiotics interact antagonistically with wildtype bacteria but synergistically with mutants (λ₁ = −1, λ₂ = 1). Here, for the sake of simplicity, we assume that antibiotics M and N have the same maximum kill rate (E_max) on wildtype bacteria (see Table 1). However, the maximum kill rate of both antibiotics declined against mutants (see Table 1) even though both same. We also assume that the IC₅₀’s of $M$ and $N$ antibiotics are the same. Nevertheless, this can be achieved by a simple change of variables, scaling by an appropirate value.

Figure 1 shows that the system (5) solution converges to the infection-free equilibrium P₀, as indicated by R_s < 1 and R_r < 1 in all scenarios (Figs. 1A–1C). Synergistic antibiotic interaction results in lower reproductive numbers for both wildtype and mutant bacteria compared to additive interaction. Conversely, antagonistic interaction leads to higher reproductive numbers. In Figs. 1D–1F, where R_s < R_r and R_r > 1, the solutions approach the equilibrium point P₁, indicating mutant bacteria evasion. In Figs. 1G–1I, where R_s > R_r and R_s > 1, the system (5) converges to the equilibrium point P₂, showing stabilization of less fit mutants by mutations from wildtype bacteria.

Impact of combination therapy on the minimum inhibitory concentration of mutants

Here we inspect the relationship between the interaction parameters of antibiotics for wild-type bacteria (λ₁) and mutants (λ₂) on the MIC. First, we analytically investigate the influence of (λ₂) on the MIC of the mutants. Therefore, we take R_r = 1, for no visible growth, in Eq. (8):

(25) $R_r={\displaystyle \frac{{\beta }_r}{\left({\alpha }_{21}{\alpha }_{22}{\lambda }_2+{\alpha }_{21}+{\alpha }_{22}\right)+{\mu }_r}=1.}$

Solving for (λ₂) we get:

(26) ${\lambda }_2={\displaystyle \frac{{\beta }_r-{\mu }_r-{\alpha }_{21}-{\alpha }_{22}}{{\alpha }_{21}{\alpha }_{22}}={\displaystyle \frac{{\beta }_r-{\mu }_r-\left({\displaystyle \frac{E_{max}^r}{I{C}_{50}^{M,R}}{\displaystyle \frac{{\theta }_1}{{\mu }_1}}}\right)-\left({\displaystyle \frac{E_{max}^r}{I{C}_{50}^{N,R}}{\displaystyle \frac{{\theta }_2}{{\mu }_2}}}\right)}{\left({\displaystyle \frac{E_{max}^r}{I{C}_{50}^{M,R}}{\displaystyle \frac{{\theta }_1}{{\mu }_1}}}\right)\left({\displaystyle \frac{E_{max}^r}{I{C}_{50}^{N,R}}{\displaystyle \frac{{\theta }_2}{{\mu }_2}}}\right)}.}}$

For the simplicity, we assume that both antibiotics have the same minimal inhibitory concentration ( IC^M,R₅₀ = IC^N,R₅₀ = IC^R₅₀), which can also be obtained by changing the variables and maximum killing rate (E^M,R_max = E^N,R_max = E^R_max) for mutants. Here we investigate only the condition where μ₁ = μ₂ and θ₁ = θ₂. Then, Eq. (26) is written as:

(27) ${\lambda }_2={\displaystyle \frac{{\beta }_r-{\mu }_r-2\left({\displaystyle \frac{E_{max}^r}{I{C}_{50}^R}{\displaystyle \frac{{\theta }_1}{{\mu }_1}}}\right)}{{\left({\displaystyle \frac{E_{max}^r}{I{C}_{50}^R}{\displaystyle \frac{{\theta }_1}{{\mu }_1}}}\right)}^2}.}$

We have found (see Supplemental Information equation S26–S30):

(28) $$I{C}_{50}^R={\displaystyle \frac{E_{max}^{r}MI{C}_r}{{\beta }_r-{\mu }_r}}$$where MIC_r is the MIC for mutants treated with a single antibiotic. Substituting Eqs. (28) to (27), we reach

(29) ${\lambda }_2={\displaystyle \frac{MI{C_r}^2{{\theta }_1}^2-2MI{C}_r{\mu }_1{\theta }_1}{\left(({\beta }_r-{\mu }_r\right){\mu }_1{)}^2}.}$

Equation (29) establishes a connection between λ₂ and MIC_r, indicating that λ₂ depends on the minimum effectiveness of antibiotics when used individually, as reflected by the maximum MIC value. This behavior is illustrated in Fig. 2, considering the values provided in Table 1. The figure demonstrates that as the interactions shift from antagonistic to synergistic, the antibiotics can have higher MIC values. Synergistic interactions compensate for the inefficiency of a single antibiotic, whereas antagonistic interactions require the antibiotics to be highly effective when used individually.

Linking antibiotic interaction to the growth rate of the mutants which capture resistance

Antagonistic interactions benefit both wildtype bacteria and mutants, prompting us to explore their evolutionary implications. By comparing the growth rates of wildtype and mutant bacteria, we can determine the dominant population (Kotil & Vetsigian, 2018). Those with a growth advantage can exert control over the gene pool, rapidly passing on their advantageous qualities to future generations. This enables the population to expedite its development when beneficial traits, such as resistance, spread and the bacteria adapt to capitalize on this advantage.

Now we suppose a quasi-stable condition for concentrations, and we compute the maximum growth rate of wildtype and mutant bacteria when s and r are close to zero in system (5). To calculate the growth rate of wildtype and mutant bacteria, we divide the Eqs. (5a) and (5b) with s and r, respectively. We reach

(30a) $G_s={\beta }_s-m-\left(\left({\alpha }_{11}{\alpha }_{12}{\lambda }_1+{\alpha }_{11}+{\alpha }_{12}\right)+{\mu }_s\right)$

(30b) $G_r={\beta }_r-\left(({\alpha }_{21}{\alpha }_{22}{\lambda }_2+{\alpha }_{21}+{\alpha }_{22}\right)+{\mu }_r)$where G_s and G_r are the growth rates of wildtype and mutant bacteria in quasi-stable conditions, respectively. To simplify the analysis, we assume that both antibiotics have the same minimal inhibitory concentration for mutant bacteria (IC^M,R ₅₀ = IC^N,R ₅₀ = IC^R₅₀) and wildtype bacteria (IC^M,S ₅₀ = IC^N,S ₅₀ = IC^S ₅₀). We also assume that both antibiotics have the same maximum kill rate for mutant (E^M,R _max = E^N,R _max = E^R _max) and wildtype (E^M,S_max = E^N,S_max = E^S_max) bacteria. So that

${\alpha }_{11}={\displaystyle \frac{E_{max}^S}{I{C}_{50}^S}{\displaystyle \frac{{\mu }_1}{{\theta }_1},{\alpha }_{12}={\displaystyle \frac{E_{max}^S}{I{C}_{50}^S}{\displaystyle \frac{{\mu }_2}{{\theta }_2}}}}}$and

${\alpha }_{21}={\displaystyle \frac{E_{max}^R}{I{C}_{50}^R}\frac{{\mu }_1}{{\theta }_1},{\alpha }_{22}=\frac{E_{max}^R}{I{C}_{50}^R}\frac{{\mu }_2}{{\theta }_2}.}$

Here, we will investigate only the condition where μ₁ = μ₂ and θ₁ = θ₂. Consequently, α₁₁ = α₁₂ and α₂₁ = α₂₂ then Eqs. (30a) and (30b) has become

(31a) $G_s={\beta }_s-m-\left(\left({{\alpha }_{11}}^2{\lambda }_1+2{\alpha }_{11}\right)+{\mu }_s\right)$

(31b) $G_r={\beta }_r-\left(({{\alpha }_{21}}^2{\lambda }_2+2{\alpha }_{21}\right)+{\mu }_r).$

Since antibiotics have less effect on resistant mutant bacteria than wildtype, so we can write α₁₁ = δ α₁₂ for some δϵR, substituting this in the Eq. (31b) we get:

(32) $$G_r={\beta }_r-\left({\delta }^2{{\alpha }_{11}}^2{\lambda }_2+2\delta {\alpha }_{11}\right)+{\mu }_r).$$

Solving the Eq. (31a) for α₁₁ we find

(33) $${\alpha }_{11}={\displaystyle \frac{-1+\sqrt{-{\lambda }_1{\mu }_s+{\lambda }_1{\beta }_s-{\lambda }_1{G}_s-{\lambda }_{1}m+1}}{{\lambda }_1}}$$substituting the Eq. (33) into the Eq. (32)

(34) $$\begin{array}{l}{G}_r={\beta }_r-\frac{{\delta }^2{\left(-1+\sqrt{-{\lambda }_1{\mu }_s+{\lambda }_1{\beta }_s-{\lambda }_1{G}_s-{\lambda }_{1}m+1}\right)}^2{\lambda }_2}{{\lambda }_1{\lambda }^2}\\ -\frac{2\delta \left(-1+\sqrt{-{\lambda }_1{\mu }_s+{\lambda }_1{\beta }_s-{\lambda }_1{G}_s-{\lambda }_{1}m+1}\right)}{{\lambda }_1}-{\mu }_r).\end{array}$$

In Fig. 3, we analyze the growth rate surface (G_r) of resistance-acquiring mutants by plotting it against antibiotic interaction parameters (λ₁ and λ₂) for wildtype and mutant bacteria. The parameter values from Table 1 are used, and G_s is set to 0 to represent complete inhibition of wildtype bacteria. Figure 3 confirms our expectations by demonstrating that G_r decreases with increasing λ₁ antagonism in wildtype bacteria and increases with increasing λ₂ antagonism in mutant bacteria. Notably, G_r is more influenced by λ₁ than λ₂. Additionally, a numerical simulation in Fig. 4 validates our analytical findings from Fig. 3.

In Fig. 4A, we simulate the temporal progression of wildtype and mutant strains under four antibiotic combination therapies. Firstly, antibiotics M and N exhibit antagonistic effects on the wild type and synergistic effects on the mutant. Secondly, their combination impacts the mutant antagonistically and the wild type synergistically. Thirdly, they have antagonistic effects on both mutant and wild-type cells. Ultimately, antibiotics M and N demonstrate synergistic effects on both cell types.

Figure 4A reveals that mutant bacteria reach their highest population density at λ₁ = 1 (synergistic). Despite the antagonistic interactions between antibiotic pairs for mutant bacteria (λ₂ = −1), the mutant population acquiring resistance is lower when antibiotic concentrations are low in the second and third combination therapies. As the antibiotic concentration approaches its maximum, the population decreases in the same combination therapy scenario. Consequently, Fig. 4A indicates that λ₂ has a minimal impact on the growth rate of mutant bacteria acquiring resistance.