THE CONTINUOUS MORBIDOSTAT: A CHEMOSTAT WITH CONTROLLED DRUG APPLICATION TO SELECT FOR DRUG RESISTANCE MUTANTS

. The morbidostat is a bacteria culture device that progressively increases antibiotic drug concentration and maintains a constant challenge for study of evolutionary pathway. The operation of a morbidostat under serial transfer has been analyzed previously. In this work, the global dynamics for the operation of a morbidostat under continuous dilution is analyzed. The device switches between drug on and drug oﬀ modes according to a simple threshold algorithm. We prove the extinction and uniform persistence of all species with both forward and backward mutations. Numerical simulations for the case of logistic growth and the Hill function for drug inhibition are also presented.


1.
Introduction. Antibiotic drug resistance is a global health problem [15]. Today, clinically important bacteria are characterized by their resistance to single or multiple drugs. Historically, penicillin-resistant Staphylococcus aureus was discovered soon after the introduction of penicillin in clinical environments [3] and still up to now the antibiotic drug resistance is still a subject of intense research [2,7,8,9,17,23]. Although the in vivo evolution of drug resistance has been proven to be useful [2], these experiments are done retrospectively and control over the environmental factors is elusive. Recently, a more advanced chemostat [20] known as a morbidostat has been developed with the goal of imposing well controlled temporal drug profile to study the evolutionary pathway [13,16,18,20]. In short, a morbidostat is a microbial selection device that progressively increase the antibiotic concentration to continuously challenge the microbe. During the course of evolution, drug resistant mutants prevail because of antibiotic drug resistance. In the end, daily frozen samples are analyzed retrospectively to reveal the genetic mutation events that lead to antibiotic drug resistance. In one incarnation demonstrated Figure 1. Schematic of a continuous morbidostat. There is no drug injection when the total microbes are less than threshold U . There is continuous drug injection once the total microbes reach the threshold U .
by Kishony's group [20], the microbial population is diluted periodically and this situation has been analyzed mathematically. In another incarnation demonstrated by Yang's group [16], the microbial population is diluted continuously and a simple threshold algorithm is used and this situation is analyzed in this work. Figure 1 shows the schematic of a continuous morbidostat device.
In [3], the cultivation of the microbes is assumed to be under serial transfer, and a simple threshold algorithm is used as feedback. In this work, we studied the case with continues dilutions and injections instead of periodical dilutions and injections. More preciously, there is no drug injection when the total microbes are less than a threshold U . It implies that without drug injection, the drug P simply decays exponentially, which is dP dt = −DP.
Thus, the drug P will die out eventually after sufficiently long time without drug injection. However, we assume that once the total microbes reach the threshold U , there is continuous drug injection in the model, which is dP dt = −D(P − P 0 ).
Hence the drug concentration P will converge to the input drug concentration P 0 under continuous drug injection. The operation of morbidostat can therefore be thought as switching between these two modes and the feedback threshold algorithm decide when the switching happens.
2. Description of our models. In the simplest scenario, we could formulate the transition from a wild type bacteria population to N mutant strains. For the drug on-drug off model with only forward mutations (see Figure 2) and both forward and backward mutations (see Figure 3), the growth dynamics with the nutrient substrate S under the influence of the drug inhibitor P are given by (1) and (2), where i = 1, 2, ..., N − 1. v 0 = u and v i are the volume densities of the wildtype and mutant populations, respectively. γ denotes the yield constant, which reflects the conversion of nutrient to bacteria. In Equations (1) and (2), g 0 (S) and g i (S) denote the growth rates of the wild type u and mutants v i , respectively. Furthermore, the growth rates satisfy g 0 (0) = 0, g 0 (S) > 0, g i (0) = 0, and g i (S) > 0 for i = 1, 2, · · · , N , which implies the fact that the bacteria do not grow if there is no nutrient substrate in the bioreactor, and a higher concentration of nutrient leads to higher growth rates of bacteria. There are two important cases in the microbiology, namely, g i (S) = m i (S) and g i (S) = miS ai+S , i = 0, 1, 2, ..., N. The wildtype and mutants are assumed to consume drug while drug inhibits the growth of the bacteria. The consumption rates on drug P by the wild type bacteria u and mutants v i are denoted by h 0 (P ) and h i (P ), respectively. The consumption functions h 0 (P ) and h i (P ) are nonnegative functions and are increasing in P . The drug inhibitions for the wild-type and mutants are described by f 0 (P ) and f i (P ), which satisfy the convention that when P = 0, f 0 (0) = 1 and f i (0) = 1. Meanwhile, f i (P ) < 0 for i = 1, 2, · · · , N since drug with a higher-concentration leads to a stronger inhibition to both wild type and mutants. Considering the fact that mutants always have stronger resistances to the inhibitor than wild-type, and mutant v i always have a  Mutant v i mutates to mutant v i+1 with a forward mutation rate q ii+1 , while mutant v i+1 mutates to mutant v i with a backward mutation ratẽ q ii+1 . We have v 0 = u and i = 0, 1, 2, · · · , N − 1. stronger resistance to the inhibitor than mutant v i−1 , we have f i−1 (P ) ≤ f i (P ) for i = 1, 2, · · · , N. In other words, we have It has shown [14] that in some cases, mutants remained abundant in the population or even have greater growth rate despite the absence of inhibitor, while in others the mutant has smaller growth rate than wild-type bacteria, as expected. In this work, we study the case that mutant have same or greater growth rate compared with the growth rate of wild-type bacteria in a drug free environment. Then the wild type grows more slowly than all the mutants, and mutant v i−1 grows more slowly than the mutant v i in the inhibitor environment for i = 1, 2, · · · , N , which is our basic assumption as follows, q ik denotes the forward mutation rate from mutant v i to mutant v k , andq ik denotes the backward mutation rate from mutant v k to mutant v i . We assume that the mutation rates q ik andq ki are quite small compared with the difference of growth rates g i (S)f i (P )−g i−1 (S)f i−1 (P ) for all i = 1, 2, · · · , N , i.e., q ik ,q ki g i (S)f i (P )− g i−1 (S)f i−1 (P ) for i, k = 1, 2, · · · , N .
) be the solution of the I.V.P (3) and (4). Then the solution exists and is unique for all t ≥ 0. Furthermore, Proof. We note that although the differential equation for P is discontinuous, the solution (S(t), u(t), v 1 (t), ..., v N (t), P (t)) is continuous. We claim that the maximal interval of existence is [0, ∞). If not, then the solution will blow up in finite time. However, the solution has a priori bound, namely, P (t) ≤ P (0) .
} as long as the solution exists. By a theorem of global existence [12], the maximal interval of existence of the solution (S(t), Adding the first N + 2 equations of models (3) and (4) together respectively, and define then the following single equation is obtained, It follows at once that, , which leads to C(t) → S 0 as t → ∞. That completes our proof.
Proof. For drug on-drug off models (3) and (4), we have From Lemma 3.1, for any ε > 0, there exists large enough t 0 = t 0 (ε) > 0 such that t ≥ t 0 implies that Therefore, for t ∈ [t 0 , +∞], we have Let S * (t) be the solution of the following system, Then we have We have η is always positive from (H2). Since From Lemmas 3.1 and 4.1, the forward mutations model (3) can be reduced to the following limiting system [19,11].
We note from [19] Appendix F, [22] Theorem 1.2.2 and [10], it can be shown that the solution of original systems (3) and (4) has the same asymptotic behabiour as those in the limiting system provided the acyclic condition is satisfied and the solution of the corresponding limiting system converges to an equilibrium. We denote (S(t), v N (t), P (t)) as the solution of the limiting system (8) with the initial values S(0) > 0, v N (0) > 0, P (0) > 0. The long-term dynamics of the limiting system is our concern. We will analyse the global dynamics of the limiting system (8) from the following three cases. 4.1. Case 1: Drug off case of the limiting system. If v N (t) < U for all t > 0, then drug P follows dP dt = −DP − h N (P )v N all the time, which is the drug off case. In this case, the limiting system (8) becomes the following system, It is easy to check that P (t) → 0 in system (9) as t → ∞. Then system (9) can be reduced to the following limiting system since lim Straightforward calculation shows that model (10) has at most two equilibria, a boundary equilibrium E 1 (S 0 , 0) and a positive equilibrium , which always holds true since S * ∈ (0, S 0 ) and g N (S) > 0. From [1,11], the following result is obtained.

4.2.
Case 2: Drug on case of the limiting system. If v N (t) > U for all t > 0, then drug P follows dP dt = (P 0 − P )D − h N (P )v N all the time, which is the drug on case. In this case, the limiting system (8) could be reduced into the following limiting system, Since S + v N → S 0 as t → ∞ by lemmas 3.1 and 4.1, then system (11) can be reduced to the following limiting system, We will first prove that every periodic solution of the limiting system (12) is orbitally stable. Let if we assume (v N (t), P (t)) is a periodic solution of model (12) with period T , and (v N ,P ) is the equilibrium of the limiting system (12), then we have It is easy to check that is locally asymptotically stable and every periodic solution (if it exists) is orbitally stable, then from [4], there is no periodic solution and thus by Poincare Bendixson Theorem, we have (v N (t), P (t)) → (v N ,P ) as t → ∞. Moreover, (v N ,P ) is the unique equilibrium of the limiting system (12), which satisfies Therefore, for limiting system (11), we have the following result.
It can be checked that v * N >v N since f N (P ) < 1 and g N (S) > 0.

Case 3: oscillation of the limiting systems.
If v N (t) oscillates around U , then the drug P follows the rule dP dt = −DP − h N (P )v N if v N (t) < U at some time t, otherwise it follows the rule dP dt = (P 0 − P )D − h N (P )v N . Therefore, the inhibitor concentration oscillates in this case as the system is trying to maintain constant bacteria density v N through feedback. We also demonstrate it by numerical simulation in Section 6.
Based on the analysis in the above 3 cases, we state our main result about the drug on-drug off model with only forward mutations in the following theorem.
The biological interpretation is that the concentrations of the nutrient S and mutant v N will stay at S 0 − v * N and v * N , respectively after a long-term. While drug P will dilute out when time t is sufficiently large.
is globally attracting, in other words, we have lim t→∞ (S(t), v 1 (t), v 2 (t), · · · , v N (t), P (t)) = (S 0 −v N , 0, 0, · · · ,v N ,P ). It implies that the wild type u and mutants v 1 , v 2 , · · · , v N −1 will die out eventually, while the concentrations of the nutrient S, mutant v N and drug P will stay at It is easy to verify that v * * N <v N < v * N .
Proof. (i) It suffices to show v N (t) < U for t large. From the equation of v N and f (P ) ≤ 1 for P > 0, we have that for ε > 0 small and t large, (ii) It suffices to show v N (t) > U for t large. Since 0 < P < P 0 for t large, which is a contradiction. Similarly, if v N (t) > U for t large, then from Theorem 4.4, v N (t) →v N as t → ∞. This leads tov N ≥ U, which is also a contradiction. Hence, v N (t) oscillates around U .

Remark 1.
We conjecture that if v * * N < U <v N , the conclusion of (ii) holds. We will verify this conjecture in Section 6.

5.
Dynamics of drug on-drug off model with forward-backward mutations.
Lemma 5.1. For model (4), if v j (t) → 0 as t → ∞ for some j, 0 ≤ j ≤ N, then v k (t) → 0 as t → ∞ for all k = j, and 0 ≤ k ≤ N. This result indicates that wild type and all the mutants will go extinction if any species of them goes extinction in the long-term.
Proof. We may assume 0 < j < N, the proofs for the cases with j = 0 or j = N are the same as that for 0 < j < N. Since v j (t) → 0 as t → ∞, then for any > 0, there exists N j > 0 such that 0 < v j (t) < for t ≥ N j .
From model (4), we have It follows that It is easy to verify that from the equations of (3) and Lemma 3.1, we have |S | ≤ D + max δ≤S≤S 0 ,0≤P <P 0 {f N (P )g N (S)} S 0 , and |v j | ≤ CS 0 , where C is a constant satisfies Similarly, |v k | is bounded for k = 0, 1, 2, ..., N . Since f j (P ) < 0, and |g j (S)f j (P )| ≤ max δ≤S≤S 0 ,0≤P <P 0 {f N (P )g N (S)}, we have | d 2 vj dt 2 | ≤ M j for some M j > 0, which implies that |v j | is also bounded. Therefore, from [5], we have v j → 0 as t → ∞, since v j → 0 as t → ∞, and v j is bounded. It implies (13). Therefore, v k → 0 as t → ∞ for k = j. That completes our proof.
For the drug on drug off model with both forward and backward mutations, we have the following result.
Theorem 5.2. Consider the drug on drug off model (4), we have two cases for the long-term dynamics of the model, (i) When g 0 (S 0 ) ≤ D, the wild type bacteria u and mutants v i , i = 1, 2, 3, · · · , N go extinction, there is no inhibitor P exist in the device. The nutrient S will persist in a concentration of S 0 in the long-term. (ii) When g 0 (S 0 ) > D, the wild type bacteria u and mutants v i , i = 1, 2, 3, · · · , N coexist, and the most resistant microbe v N dominates the rest of the other species provided the mutation ratesq i are sufficiently small.
Proof. From Lemma 5.1, either all species go extinction or all species coexist.
(i) If all species go extinction, then the system (4) can be reduced to the following limiting system, Hence we have S → S 0 , and P → 0 as t → ∞. From the second equation in the drug on and off model (4), it follows that g 0 (S 0 ) ≤ D. Therefore, for the drug on drug off model (4), when g 0 (S 0 ) ≤ D, we have the wild type bacteria u and mutants v i , i = 1, 2, 3, · · · , N go extinction. Furthermore, there is no inhibitor P exist while the nutrient S will persist at the concentration of S 0 in the long-term.
6. Numerical simulations. In this section, we carry out numerical simulations with realistic parameters to verify our theoritical results.
For simplicity, we assume that both the wild type and the mutants have equal uptake function and the growth takes the logistic form, namely, g 0 (S) = g i (S) = mS, i = 1, 2, 3, · · · , N.
We assume the functions f 0 (P ) and f i (P ) take the Hill function form of order L, which are where i = 1, 2, 3, · · · , N . Note that the order L stems from the allosteric cooperativity of the drug inhibition [21]. The drug inhibition effects depend on the detailed mechanism. For example, they can result from the binding of the antibiotic drug to the metabolic enzyme, which synthesizes the key precursor of biomass production of the bacteria. Taking trimethoprim (TMP) as a specific example, this antibiotic binds to dihydrofolate reductase (DHFR), an enzyme that controls the biosynthesis of folic acid. The mutation of the gene encoding DHFR will modify the binding affinity of TMP [20]. The parameters K 0 and K i in (15) can actually be extracted from the experimental values of IC 50 , defined as the drug inhibitor concentration at which the growth rate is 50% of the maximal growth rate. The sample volume in the culture vessel of the morbidostat is 10 ml, and the confluent density of EColi is typically 10 9 cell/ml. For simplicity, we assume that the morbidostat operates at around 10% of the confluent density, which is set by the threshold population U . We can conveniently set the yield constant γ to be 1, so that one unit of substrate density will transform to one unit of bacteria. With the volume units set to 0.1 nl, the input substrate S 0 = 10 corresponds to 10 8 cell/ml.
The constants K 0 and K i in the drug inhibition function are set to 1 and 10 i , and the order in the Hill function L = 1, 2, 3. We first consider the forward mutation rates q = 10 −2 hr −1 . Typically, the dilution ratio D = 0.9. The growth rate is set by the constant m to be 0.3hr −1 in the logistic growth function. The initial conditions are S(0) = 10, u(0) = 0.01, v(0) = 0, and P (0) = 0. We demonstrate our theoretical results by conducting the following simulations.
6.1. Three mutants with only forward mutations. We first assume that there is wild type u and three mutants v 1 , v 2 and v 3 in the bioreactor and there is no backward mutation. Then we have the drug on-drug off model in the following form, Straightforward calculation shows that v * 3 = 7.000 andv 3 = 6.327 by the parameters given above and P 0 = 10, r = 0.5 and a = 0.5. We verified the three cases of our theoretical analysis in Section 4 using numerical simulations as follows.
(iii) If U = 6.5, then we havev 3 < U < v * 3 , which is Case 3 according to Section 4. From the simulation results, we can see the inhibitor concentration P oscillates between 0 andP , which is P ∈ (0, 6.728), while the bacteria density v 3 will be maintained at nearly constant around U , please see Figure 6.