Exploiting ecology in drug pulse sequences in favour of population reduction

A deterministic population dynamics model involving birth and death for a two-species system, comprising a wild-type and more resistant species competing via logistic growth, is subjected to two distinct stress environments designed to mimic those that would typically be induced by temporal variation in the concentration of a drug (antibiotic or chemotherapeutic) as it permeates through the population and is progressively degraded. Different treatment regimes, involving single or periodical doses, are evaluated in terms of the minimal population size (a measure of the extinction probability), and the population composition (a measure of the selection pressure for resistance or tolerance during the treatment). We show that there exist timescales over which the low-stress regime is as effective as the high-stress regime, due to the competition between the two species. For multiple periodic treatments, competition can ensure that the minimal population size is attained during the first pulse when the high-stress regime is short, which implies that a single short pulse can be more effective than a more protracted regime. Our results suggest that when the duration of the high-stress environment is restricted, a treatment with one or multiple shorter pulses can produce better outcomes than a single long treatment. If ecological competition is to be exploited for treatments, it is crucial to determine these timescales, and estimate for the minimal population threshold that suffices for extinction. These parameters can be quantified by experiment.

S1 Text, respectively. We see that for smaller τ , the region where the higher N -sequences get better increases. This occurs especially for negative s, where the region where the final pulses of the higher N -sequences are best (bottom right corner, marked with black line in Fig. 1 of the main text) extends towards much smaller t r , and merges with the region of the bottom left, an intermediate pulse in the higher N -pulse sequence is best. The extension towards smaller t r of the region where the last pulse is best can be understood by an argument similar to the main text discussed around Fig. 4 d-e, where one gauges what t (f ) w is required for the first pulse to be better; since τ is smaller, the t r where t r = τ − t (f ) w ≡ τ − t final o is smaller too. For larger τ , the opposite is the case for the bottom right corner, for exactly the same reason. The regions where the first pulse of longer sequences are best moves down towards more negative s, but also shrinks in t r compared to τ = 60, as for τ = 90 the t (1) r of the single pulse is here already much longer than t o + t (N ) r (see the main text). However, the regions do not change qualitatively. This might of course be different if one considers τ < t o or τ < t final o , but this would go beyond the scope of this work.  Fig. 1, for different τ : τ = 30 for the left (a) and τ = 90 for the right panel (b), respectively. We see that the regions where the pulse sequences are better shift consistently: the region where the higher pulses of the sequences yield the lowest global minima (bottom right) shifts towards the left for increasing τ , while the region where the first pulse of the sequence is better than the single pulse (top left) moves inwards and out. In general, for larger τ , the single pulse does better for a larger area in (s, t r ) space. Fig. 1 of the main text for different cost for the resistant species, switching rates and death rates

2.
In the main text, we chose a fitness cost of k = 0.1 for the growth rate of the more resistant species, since this parameter gave a good, yet not too confusing overview over where which pulse sequence yielded the best results. In Fig. B, we show the analogous figures to Fig. 1 in the main text for two different costs k = 0.4 and k = 0.01 (top panels, (a) and (b)). These different costs change the values for the optimal initial and final times spent in the low environment, t o , and t final o (see section 3 of the supporting information).
The general features of Fig. 1 remain the same: the single pulse is best for the largest part of the phase diagram and is lowest at the thick blue line, where the skewness is such that the initial time in the low environment matches the t o for this parameter combination. The first pulse of a longer pulse sequence is best in the top left corner, around the corresponding thick line where the initial time in the low environment for this first pulse matches t o . The region marked with the black curve corresponds to the region where the last pulse of these higher order pulse sequences is best. Compared to Fig. 1 in the main text, the increased cost of k = 0.4 (Fig. Ba in S1 Text) leads to an increased importance of this last region. In addition, the regions where intermediate pulses of the longer pulse sequences are best (marked with white lines), are also larger than in Fig. 1 of the main text, but occur at exactly the predicted positions (in the bottom right corner, and between the regions where first pulses of longer sequences yield lowest minima, in the top left).
For the smaller cost of k = 0.01 (Fig. Bb) compared to k = 0.1 in Fig. 1 in the main text, the behaviour is analogous to what we have described so far: the regions encircled in white disappear, and region where the last pulse of the longer sequences yield lower minima (marked with the black line) becomes smaller. In addition, similarly to what happens for longer τ (see Fig. F in S1 Text), the importance of the region where the single pulse yields best results increases. This is due to the fact that the optimal time spent initially in the low environment t o is shorter for shorter costs, and that thus τ /t o increases. Thus means that the t r of the single pulse is longer, and can reduce the  Fig. 1, for different rates and costs for mutation: a) k = 0.4, δ = 0.1; b) k = 0.01, δ = 0.1; c) k = 0.4, δ = 0.3; (so far for all µ w = 10 −6 , µ r = 0, like in the main text) d) µ w = 10 −7 , µ r = 10 −4 and δ = 0.1, k = 0.1 like in the main text. The general features are independent of the system parameters. For the largest part of skewnesses and time in the high environment, including the region where the skewness is such that initial time in the low environment matches the optimal initial time spent in the low environment, t o , for the corresponding parameters (thick blue lines), the single pulse is best. The first pulse of a higher order pulse sequence is best when the initial time in the low environment of this first pulse matches t o , and in areas around it (top left corner). The last pulse of such higher sequences is best for large t r and negative skewnesses (marked by a black line). The dominant differences between the different figures are the exact boundaries of these regions, and the regions where intermediate pulses of the longer pulse sequences yield the lowest minimum (white lines). Due to the high decrease in population for the high death rate of δ = 0.4, we could not precisely resolve the exact position of the black boundary and thus marked it with dashes. population for longer, than t o + t r /N for the N th pulse sequence.
In addition, the changed costs influence the absolute value of the minimal population at the best pulse. The minimal population then is smaller for higher fitness costs, as the population decays more strongly due to the slower growth of the more resistant species.
Finally, we also investigate a combination of increased death rate and cost (Fig. Bc, to be compared with the Fig. Ba above for the same cost), as well as different switching rates. Again, the features of the figures are exactly the same as in Fig. 1 of the main text. The increased death rate meant that the population reduction was higher, and so for high t r the different pulse sequences gave very small and effectively indistinguishable lowest minima. Thus, we could not resolve the region where the last pulse of the longer sequences was best precisely, and marked it with a dashed black line in these cases.
Changing the switching rates ( Fig. Bd) also leads to only very small changes of general features compared to Fig. 1 in the main text. We note the increase in the areas marked with the white lines, where intermediate pulses of higher pulse sequences do best, but again at positions where they would have occurred for Fig. 1. Thus, we conclude that the main features of, as well as the relevant explanations for Fig. 1 in the main text, are robust against parameter variations. Panel (a) of Fig. C shows the variation of t o as a function of death rate δ and cost k in the growth rate of the resistant species, k = λ − λ r . Obviously, t o increases when this cost increases, as the resistant species will take longer to grow. For fixed cost, there is a minimal t o at a certain death rate, with t o increasing towards both small and large death rates. This is due to the fact that for small growth rate costs for the resistant species, the death rate dominates for the wild type which decays faster, while for large growth rate costs, the resistant species' growth is also drastically limited by the high death rate. This behaviour is most pronounced for large cost k = 0.5, where t o diverges for δ = 0.5, since then r cannot grow any more at all. In general, we see that t o can acquire high high values of around t o ≈ 50. This means that the optimal t (i) w can have a very large effect, and explain why sometimes more moderate treatments are functional, even though they come with the risk of evolving resistance in principle.

Optimal t
In Fig. C b we show how t o depends on the switching rates µ r and µ w . Even if those rates do not change the actual dynamics very much, they have a big influence on t o , as they crucially influence the initial condition, that is the attractive fixed point of the (free) environment. For larger µ w the population starts out with a relatively higher amount of resistant species and less wild-type so that the resistant species can take over the population faster, thus decreasing t o . We note that t o decreases linearly with µ w for µ w up to µ w ≈ 10 −2 . The dependence on µ r is less pronounced for small µ w and µ r as to lowest order in µ w,r the fixed point of the free-system only depends on µ w .
Concerning t final o , we note that it depends on the preceding t

Comparing minima within a pulse sequence
There are two competing concepts that determine during which pulse the population minimum could occur: First, the time during which the population effectively decays exponentially, which can be larger than t r for pulses early in the sequence where w is still present, and will eventually be t r for pulses late in the sequence; and second, the initial population, from which the decay leading to the lowest minimum starts. The latter might be lowest in intermediate pulses, where w has decayed, but r has not yet taken over the population.
We want to discuss at which t r a second pulse can yield a lower minimum than a first pulse for s = 1 and s = −1.
For s = 1, we have seen that the maximal period during which the population decays is t w + t r is small, there could be another contribution from t (i),2ndpulse w , that is a contribution like t (f ) w only at the beginning of the 2nd pulse], which is maximal when t (i) w = t o . Since pulses are periodic, and since we consider τ < 126 (which effectively is a requirement on t r being such that the population will be lower than the minimum of the (low) regime (obtained after t o ) after the second t (i) w regime), n will grow up to r low in the following (low) regime of the second pulse; there, the population can only decay during t r , and thus the second pulse will have a higher minimum for t (i) w = t o . More precisely, for pulses with t (i) w > t o , the period of effective population decay in the first pulse shortens, and will correspond to max(t o , t r ) for pulses with t (i) w > t o + t regrowth , where t regrowth in our case corresponds to approximately 10. In the (low) regime of the second pulse, the population will definitely increase if t (i) w > t o , meaning that in the second pulse, the effective decay period will only be t r . The decay period of the second pulse is thus either equally long as during the first pulse (t r ), or shorter when t r < t o , with the starting population of decay being definitely at r low . In the latter case, the second pulse will definitely have a higher minimum; in the first, where t (i) w > t o + t regrowth and t r > t o and the decay periods of both pulses are effectively the same, the second pulse will have an identical minimum to the first pulse. (Also, this latter case only applies for sequences with τ ( N ) > 2 * t o + t regrowth ≈ 35). Thus, the second pulse has a higher or equal minimum whenever t (i) w > t o . For s = −1, we already noted one can have a very low minimum in the second pulse if the population decays consistently for t r . However, the second pulse will actually yield a lower minimum even for larger t (f ),(N ) w : unless t (f ),(N ) w is larger than the time r needs to grow up to its fixed point r low , the effective decay will be longer than t We can estimate at which t r we expect the n(t = t r ) = w * (iii) e −δτ + r * (iii) e −δτ e λrt off to be lower than n * (iii) : this occurs for t (f ) w < (1/λ r ) ln(n * (iii) /m * (iii) ) + τ δ/λ r , which corresponds to t off < 16 for the sequence with two pulses where 2τ (N =2) = τ = 60, and thus 2t Fig. 1 shows that for t r > 40 the sequence with two pulses performs better than the single pulse, as here the second pulse of the two pulse sequence is lower.
For this t r , a small increase in s from s = −1 implies the presence of the t (i),(N ) w , which means that r regrows faster during the t Thus, the first pulse in the sequence will yield a lower minimum again, and then certainly the single pulse is better.
5. Phase space trajectories (Fig. 5 a) for skewnesses s = −0.5, s = 0.2 and s = 0.5 In Fig. 5 a we show phase space trajectories for pulse sequences with τ = 60, t r = 10 and skewness s = 0.9, consisting of one, two or three pulses, respectively. The trajectories start at the (free) fixed point close to the w-axis and evolve towards the r-axis. r/w increases from pulse to pulse since the more resistant species r steadily takes over the population. Furthermore, the population composition at the end of the pulse sequence is similar for all three trajectories. The same holds true for other skewnesses. In  figure) and r/w steadily increases from pulse to pulse. Furthermore, the population composition at the end of the pulse sequences is similar for one, two or three pulses.

Constant amount of antibiotics: Estimate
As shown in Fig. E, keeping t sum = τ + t r constant shows that the dependency on the skewness is much less pronounced than for keeping τ constant, which we explain in the following.

Fig E.
Population minimum n min (a) and composition r/w (b) for a single pulse of fixed τ + t r instead of constant τ and fixed t r . Both quantities show little dependency on the skewness s. For this case, there are two competing effects: the smaller t r , the less the resistant species is depleted, but simultaneously the smaller t r the larger τ and thus, the wild-type is depleted longer. As a result, not only optimising the onset but also avoiding a large offset is beneficial.
For the case of fixed t sum , there are two competing effects: on the one hand, for small t r , τ = t sum − t r is larger than for large t r , and thus for small t r the population w can be depleted longer. On the other hand, larger t r increase the chance to reduce the resistant species r, and thus also n. Thus, for s = 1, where the pulse starts with the (low) regime (see Fig. 2 in the main text), it makes sense to exploit the full optimal onset time, such that τ = t r + t o , since then r has not grown much until the beginning of the (high) regime. This then means that t sum = 2τ − t o , and thus τ ≈ t sum /2 + 7.7. On the contrary, for s = −1 when the pulse starts with the (high) regime, it is best to use t w = 0, large t r ) in order to prevent the resistant species from taking over the population in the subsequent (low) regime. So, it is ideal to use τ ≈ t sum /2 + t final o (t (i) w = 0, large t r )/2 ≈ t sum /2. Together, these two strategies mean the pulses are optimised in such a way that their durations do not differ drastically for the different skewnesses. Thus, the best minimum and the mutant/wild-type ratio is approximately the same for all skewnesses with a little improvement (deterioration) towards skewness s = 1 for the minimum (ratio). This is also reflected in the ratio r/w in FigEb, which sharply increases until it makes sense to use t r > 0, namely at τ + t r = t o . For low s, the ratio r/w is slightly higher, as the (low) regime where the resistant species can grow, is at the end of the pulse. For higher values of τ + t r , the ratio is very similar across all skewnesses.