Comparing optimization criteria in antibiotic allocation protocols.

Clinicians prescribing antibiotics in a hospital context follow one of several possible 'treatment protocols'-heuristic rules designed to balance the immediate needs of patients against the long-term threat posed by the evolution of antibiotic resistance and multi-resistant bacteria. Several criteria have been proposed for assessing these protocols; unfortunately, these criteria frequently conflict with one another, each providing a different recommendation as to which treatment protocol is best. Here, we review and compare these optimization criteria. We are able to demonstrate that criteria focused primarily on slowing evolution of resistance are directly antagonistic to patient health both in the short and long term. We provide a new optimization criteria of our own, intended to more meaningfully balance the needs of the future and present. Asymptotic methods allow us to evaluate this criteria and provide insights not readily available through the numerical methods used previously in the literature. When cycling antibiotics, we find an antibiotic switching time which proves close to optimal across a wide range of modelling assumptions.

Introduction

Throughout the twentieth century, a small number of discoveries have fundamentally reshaped our world. Transistors have created a world of computation and communication [1-3], the Haber–Bosch process for nitrogen fixation has massively increased our ability to produce crops [4,5], and the development of antibiotics [6,7] has not only redefined our battle against disease but also acted as a key enabling technology in surgery and intensive care [8-10]. Antibiotics have helped lift life expectancy from 47 at the beginning of the twentieth century, to the 79-year average we see today [11].

At the same time, evolution, that same process that gave birth to penicillin, ensures that antibiotics are a limited resource; every time a given antibiotic is used to save a life or prevent an infection we inevitably select those bacteria most able to resist for survival. Over weeks and years this sustained selective pressure has given rise to a variety of multiresistant ‘superbugs’ [12]. Evolution of multiresistance is further accelerated by horizontal gene transfer (HGT) [13,14], a process by which one bacterial lineage may share genetic material with another. Bacteria have, over time, developed many forms of antibiotic resistance (ABR), first to penicillin, and then to erythromycin, methicillin, vancomycin and carbapenems [15,16].

While a large number of antibiotics have been developed over the past century, the vast majority act through only a small number of essential mechanisms; as an example amoxicillin, cephalexin, doripenem, meropenem, aztreonam, ceftolozone and many others all act through the same β-lactam group as penicillin [17]. Bacteria that are resistant to one of these compounds will quickly develop enhanced resistance to the others [18]. In the past 30 years, only a few truly novel antibiotic compounds have been developed (for example, complestatin and corbomycin, discovered 2020, [19]), and of those compounds discovered, even fewer have made it to market (see table II of Coates et al. [20] for a list of antibiotics discovered in the past decades, and where they are in the process of clinical trials).

These practical difficulties are further exacerbated by the fact that antibiotic research is unprofitable [21,22]; development of new antiobiotics is hugely costly, yet any new drug is liable to be kept on hospital shelves as a ‘reserve’ antibiotic, and if it is used, will cure a patient with days or weeks, as opposed to the years or lifetime for drugs designed to treat chronic illnesses, such as cancer, depression or high blood pressure.

Given the substantial human and economic costs caused by resistant and multiresistant bacteria in clinical settings [23], currently including billions of dollars, and tens of thousands of lives, there is substantial interest in understanding how best to slow the evolution and spread of multiresistant genotypes. While there is substantial interest in developing ‘evolution proof drugs’, targeting bacterial virulence [24], using chemical mechanisms to inhibit resistance mechanisms directly [25], or pairing antibiotics where resistance to one drug implies sensitivity to another [26], for the time being, our primary tool in minimizing ABR is simple antibiotic stewardship: selecting the time and place of antibiotic use so as to minimize selection for resistance.

Antibiotic stewardship has been studied from both a clinical point of view (see the review articles by Dik et al. [27] and Chatzopoulou et al. [28] as well as [29-31]) and from a theoretical standpoint (see the literature described below).

The archetypical question asked in this context is ‘when and how should clinicians prescribe antibiotics?’ Is a hospital better off prescribing a variety of different antibiotics or applying a monoculture, switching the drug of choice every few weeks? Should all patients be prescribed multiple antibiotics (combination therapy) so as to ensure recovery, or is a lighter touch better in the long run? And how, precisely, do we define ‘better’ anyway? Is it more important to ensure optimal patient outcomes in the present, or slow the development of multiresistant bacteria over an evolutionary timescale?

A first inquiry into these conundrums was made in a simulation study by Bonhoeffer et al. [32]. Subsequent models have built on Bonhoeffer’s approach, exploring horizontal gene transfer [33,34], extensions to a larger number of antibiotic variants [35,36] and the effects of stochasticity [37]. Ibargüen-Mondragón et al. use control theory to consider both the effects of antibiotic treatment and immune response [38], while Houy et al. consider the question from the point of view of ‘what information do decision makers have access to?’ [39]. McCloud et al. explore the importance of linkage disequilibrium, the tendency of two resistance genes to be found together or apart, and model evolution on the scale of multiple inter-connected populations (for the sake of our discussion, one might imagine patients in multiple different hospitals) [40].

Concerningly, among these (and many other) investigations, no consensus on optimal treatment protocols has been reached. Uecker & Bonhoeffer [41] explain these contradictions in their clear, concise and comprehensive review article, in which they provide an overview of the subtle modeling decisions and mismatched optimization criteria that play into various contradictory rankings of antibiotic deployment strategy.

The focus of this article will be following up on a number of questions raised in Uecker & Bonhoeffer’s review. We compare the strengths and weaknesses of the various proposed optimization criteria in the literature. We propose a new, more physically justifiable optimization criteria, incorporating both patient health and evolutionary risk, a middle path between the various criteria considered previously. In addition, we find a number of analytic results that will hopefully reduce the need for time-consuming numeric exploration of parameter space, and give somewhat clearer insight into how parameter values combine to determine patient outcomes. For the important case of cyclic treatment protocols, we use our asymptotic results to identify a ‘saturation time’ tsat which gives ‘almost optimal’ cycle time across a wide range of assumptions. We observe that optimization criteria focused on delaying multiresistance often harm patient health in both the short and long term: they are often maximized by selecting treatment protocols with patient outcomes worse than the long-term effects of antibioitic resistance.

Our goal in what follows is not to question the orthodoxy of aggresive antibiotic treatment [42]. Instead our sole goal is to examine the methods by which we evaluate antibiotic management protocols, and explain the root cause of a number of conflicting and contradictory results previously published.

In §2, we introduce the basic model and current antibiotic protocols as described in the literature. Section 2.2 introduces four possible optimization criteria, both novel and historic. In §3, we examine the mean number of uninfected patients for a variety of antibiotic protocols and gather several novel analytic approximations for patient health. In §4, we explore the behaviour of four different optimality criteria, and how ‘optimal’ results change (or don’t) depending on which criteria is used and how multiresistance arises. We demonstrate how a number of optimization criteria are directly antagonistic to one another, and indeed, to patient health. These results are summarized in §5. We find combination therapy to be appropriate over a wide range of contexts, and also identify tsat, a cycle time that is ‘close to optimal’ over a wide range of modelling assumptions. It is our hope that by identifying why different papers have reached such contradictory conclusions, we can guide medical practitioners and future modellers towards better criteria, or at the very least, give them the tools to evaluate which criteria are most relevant.

Model

Let us begin by introducing a small amount of biology and defining the two models of in-hospital ABR spread that we will be building upon. Primarily, we make use of Bonhoeffer et al.’s original model from 1997 [32]. This model has acted as the basis of many of the papers that came after, and hence acts as an excellent testing space for comparing and understanding the various optimization criteria found in the literature. The second model we considered was proposed more recently by Uecker & Bonhoeffer [43] and is designed to overcome one particular simplification in the original model, which will be relevant in our future discussions. This simplifying assumption, and when it is and is not appropriate, will be discussed later.

These models are selected primarily for their analytic tractability, clarity, and foundational position in the literature. Our goal in this paper is to compare differing optimization criteria and resolve conflicting recommendations as cleanly as possible, Bonhoeffer’s model is perfect for this. With that said, future researchers may find that more detailed models incorporating more recent clinical understanding of effects, such as horizontal gene transfer [44,45] and bystander selection [46], may be more appropriate when making treatment recommendations. With that said, let us start by describing Bonhoeffer et al.’s 1997 model [32].

Imagine a hospital or hospital ward containing a number of patients. There exists some bacterial infection (for example, Streptococcus pneumoniae or Staphylococcus aureus) that spreads through the patient population and can be treated using one of two frontline antibiotic drugs (such as linezolid or telavancin). For the sake of generality, we refer to these drugs as A and B. Bacterial infections will either be susceptible to treatment, resistant to one or other of the available antibiotics, or doubly resistant. Denote the population of patients infected with each of these classes of infection as S, R, R and R, respectively. The hospital also includes a population of uninfected patients, X, who are nonetheless ‘exposed’ to infection. These patients represent patients recovering from surgery (or similar medical conditions unrelated to infection). While infection may provide an additional load, uninfected status does not correspond to a clean bill of health, and conversely, an infected patient may still be healthy enough for discharge if they have recovered from surgery, chemotherapy, or whatever was their primary cause for entering the hospital. It is also important to note that here ‘susceptible’ and ‘exposed’ are not used in the same manner as standard epidemiological SEIR models. Here ‘susceptible’ refers to the infecting bacteria’s susceptibility to antibiotic treatment, and ‘exposed’ refers to a patient's proximity to possible infection, as opposed to the presence of bacterial infection that has not progressed to the ‘infectious’ stage, as implied in SEIR models.

With these five classes of patients, we now construct a compartment type model [47]. Patients switch from one infected status to another via a variety of processes (figure 1). Patients arrive at the hospital at some rate m, m, m, m and m; immigration rate for each class will depend on the prevalence of ABR in the community. Patients in all compartments are discharged at some rate μ; this rate is assumed to be equal across all compartments. In appendix A, we relax this simplifying assumption and consider a model that explicitly differentiates between discharge due to death or recovery. The exact details of discharge have minimal impact on all major results, and hence we assume discharge rate μ is independent of infection status, as is traditional in the literature.

Figure 1

Schematic diagram of the compartment model, as defined in equations (2.1a)–(2.1e). The diagram shows the five model compartments and the flow of individuals between them. Mutation events are marked with dashed lines and are treated as being rare stochastic events. More common infection, recovery, immigration and emigration events are indicated by continuous lines; differing line thickness is used to indicate that some events are more common than others (for example, β > β, τ + γ > γ). Line thickness is not to scale. Depending on the context, the R compartment of the above model will, or will not, be included in the analysis.

Exposed individuals become infected according to mass action kinetics at rates β S X, β R X, β R X and β R X. Differences in infection rate β represent the metabolic cost a bacteria must pay in order to maintain resistance; more resistance comes at a higher cost, hence β ≥ β, β ≥ β. All β are assumed to be fairly similar, differing by only 1 or 2%, as opposed to an order of magnitude. Hospital policies that involve isolating ABR infected patients may be used to artificially reduce β further [48]. While these strategies can be very effective, they will not be our focus here.

Infected individuals recover naturally at some rate γ, and recover at some higher rate τ + γ if administered suitable antibiotic treatment (hence, administering drug A leads to recovery rate τ + γ in the S and R population and recovery rate γ in the R and R population.) We model which antibiotics are currently being administered via the indicator functions χ(t) and χ(t). These functions take the value zero or one, depending on whether or not the antibiotic in question is being administered at time t.

Single resistance strains are assumed to be prevalent enough in the community that de novo evolution of single resistance can be treated as negligible. The validity of this assumption will depend on the population under study.

Finally, multiresistance can be introduced to a hospital either via de novo mutations, horizontal gene transfer or importation from the wider community. Each of these mechanisms leads to different predictions and recommendations, which will be discussed in more detail on §4. For the time being, we consider two separate cases: the behaviour of the system either before multiresistance (R(t) = 0 = m), and the behaviour of the system after multiresistance(R(t) > 0). The transition from R = 0 to R > 0 is discussed in §4.

Taken together, these processes lead to the following five compartment model:

The above system of equations assumes that the same antibiotic regime is applied to all patients: χ(t), χ(t) take the values zero or one. These assumptions are inappropriate however when administering different drugs to different patients. The most obvious way to represent intermediate prescription values in the above model would be to select χ(t) = 0.5 = χ(t); unfortunately, this would correspond to assigning all infectious patients both drugs half the time (for example drug A in the morning, and drug B at night). In clinical practice, this does not happen, and instead ‘mixed’ drug regimes refer to the practice of assigning half of the patients drug A across the entire course of their treatment, and the remaining patients drug B.

In order to properly model this, we use the seven compartment model of Uecker et al. [43]). In this model, we track both the resistance status of infections and the prescription status of the corresponding patients, splitting compartment R into and R into . Here, subscripts represent the resistance profile of the infection, while superscripts represent the drug currently prescribed. Hence, and represent effective treatments, while and represent ineffective treatments. Treatments are ‘corrected’ at a rate q, transferring patients from and . Because all prescriptions are assumed to be equally effective against susceptible bacteria, there is no need to split the S compartment, similarly with doubly resistant infections. χ(t) and its complement χ(t) no longer represent the probabilities of receiving a particular drug in the present, but instead the probabilities of being referred to a particular treatment group. The governing equations for and arewith similar equations governing and . and are as defined in equations (2.1a) and (2.1e). We also consider an antibiotic switching rate q; this corresponds to the rate at which ineffective antibiotics are replaced by effective antibiotics in the case of single resistance. In the limit q = 0, patients are kept on their initial prescription indefinitely no matter what. In the limit q → ∞, patients are shifted to optimal treatment immediately. Parameter q can be thought of as a proxy for the intensity of testing for ABR within a given hospital system. A schematic of this behaviour is given in figure 2.

Figure 2

Schematic diagram of Uecker’s 7-box model. In this model, we track the drug being applied to individual patients: R is split into and (A resistant bacteria being treated with A or B, respectively). Although S bacteria are treated with either one drug or the other, there is no need to track which (assuming similar recovery rates under treatment). Similarly with doubly resistant infections. Ineffectual treatment combinations ( or ) are replaced by effective treatment options ( or ) at some rate q (drug switching, turquoise arrow). Differences in line thickness are indicative of differences in the corresponding rate constants (recovery from is faster than recovery from , for example).

In what follows, Uecker’s 7-box model is used when 0 < χ < 1 and Bonhoeffer’s 5-box model is used when χ ∈ {0, 1}. See table 1 for a list of parameters, their descriptions and (where appropriate) default values.

Table 1

Description of parameters, their standard values (where appropriate) and a summary of optimization criteria used throughout the main text.

description	symbols	default value
number of patients on the ward but not currently infected	X(t)	—
patients infected with AB susceptible infections	S(t)	—
patients infected with drug A resistant strain	RA(t)	—
patients infected with drug B resistant strain	RB(t)	—
patients infected with multiresistant strain	RAB(t)	—
drug A resistant patients treated with drug A	RAA(t)	—
drug A resistant patients treated with drug B	RAB(t)	—
drug B resistant patients treated with drug A	RBA(t)	—
drug B resistant patients treated with drug B	RBB(t)	—
infection rate of susceptible strain	βS	1
infection rate of drug A resistant strain	βA	0.99
infection rate of drug B resistant strain	βB	0.99
infection rate of multiresistant strain	βAB	0.98
arrival rate of uninfected patients at hospital	mX	∼2/45 days
arrival rate of patients with non-resistant infections	mS	0.5 mX
arrival rate of drug A resistant patients at hospital	mA	0.002 mX
arrival rate of drug B resistant patients at hospital	mB	0.002 mX or 0.004 mX
arrival rate of patients with multiresistant infections	mAB	0
hospital discharge rate (μ = mX + mS + mA + mB)	μ	1/15 days
baseline recovery rate	γ	1/10 days
recovery rate with suitable treatment	γ + τ	1/2.5
rate at which inappropriate drugs are replaced by appropriate drugs (seven box model only)	q	0
cycle time for drug cycling	T	varies
proportion of patients receiving drug Aa	χA(t)	varies
proportion of patients receiving drug Ba	χB(t)	varies
Rate of mutation. Multiplied by various different things depending how RAB is assumed to arrive.	ν	1 or 0.1
total uninfected ‘person days’ in first year	X365	∫0365X(t) dt
time until multiresistant arrival	Tε	—
time until half of all infections are multiresistant	T1/2	—
total uninfected ‘person days’ until T1/2	XT	∫0T_1/2X(t) dt
long-term average X(t) value, prior to RAB arrival	X¯	∫t_0t_0+TX(t)T dt
long-term average X(t) value, after RAB arrival	X¯_AB	∫t_0t_0+TX(t)T dt
smoothed ‘gain in X’ up until Tε	XT*	(X¯−X¯_AB)T_ϵ

aInterpretation of these parameter differs between the 5-box and 7-box model.

Common antibiotic management protocols

In this article, we will compare three main antibiotic management protocols frequently discussed in the literature: combination therapy, mixing and cycling (mixing and combination therapy are most common in clinical practice). These protocols define which antibiotic is initially prescribed to a patient when they are admitted to the hospital, or after surgery. This prescription may later be changed based on either patient recovery (or lack thereof), or when ABR infection is detected (see, for example, the Dutch ‘search and destroy’ policy for ABR [49]).

The simplest of these protocols, combination therapy, prescribes both A and B to all infected patients at all times. This approach is intended both to improve patient outcomes and to prevent multiresistance from arising by eradicating single resistant strains as quickly as possible. Unfortunately, combination therapy leads to increased direct costs and potentially heavier side effects. It also increases total drug prevalence, leading to concerns that it may encourage broad spectrum antibiotic resistance [50]. This duel action of increasing total antibiotic prevalence, but rapidly quashing single resistant strains, leads to some uncertainty with respect to the net effect of combination therapy on ABR. In this study, combination therapy is represented using Bonhoeffer’s 5-box model, with χ(t) = χ(t) = 1.

Mixing protocols assume that each patient is assigned either drug A or B with some probability, usually (but not always) χ = χ = 0.5. While many papers [33,36,51] have studied mixing using Bonhoeffer’s 5-box model (or similar), Uecker’s more detailed 7-box model is more faithful to clinical reality and will be the model used here whenever mixing is discussed.

Cycling protocols treat all patients with the same drug at any point in time, and switch back and forward between two (or more [35]) treatments every T days, preferentially treating with drug A for the first T days, and with drug B for the next. These time periods are generically (though not always [52]) assumed to be equal. Mathematically, cycling is represented bywith χ(t) = 1 − χ(t). Both the 5-box and 7-box model can be used to represent cycling; which is more realistic will depend on the exact implementation of cycling used in clinical practice. For the sake of analytic accessibility, we study cycling in the context of the 5-box model. Basic simulation experiments indicate that the difference between the two models is negligible, except in the case of short cycle times, where fast cycling behaves like mixing.

Other, more detailed, management protocols have been studied. Beardmore & Peña-Miller [53] make use of detailed control theory techniques in order to determine optimal aperiodic antibiotic rotation protocols, cutting through the more heuristic approaches used elsewhere in the literature. Kouyos et al.’s [37] numeric exploration expands the space of possible treatment protocols by considering ‘informed switching’ protocols adapted for the stochastic hospital environment. Consideration of these more complex protocols is beyond the scope of this paper, but we do recommend these past works to the interested reader.

Optimization criteria

Broadly speaking, antibiotic protocols seek to achieve two conflicting goals: to maximize patient health and minimize the rate at which resistant bacteria (especially multiresistant bacteria) arise and develop. While these goals are easy enough to understand in an intuitive sense, there have nonetheless been a number of different formulations mathematically; each with their own strengths and weaknesses. For the sake of brevity, we can not consider all possible optimization criteria ever considered in the literature; our goal here is to present a representative sample. Additional optimization criteria, such as exponential discounting (as used in, for example [54]) and Bonhoeffer’s G1/2 are discussed in appendix D and B.4, respectively.

The most straightforward evaluation criterion is simply to maximize the number of ‘healthy’ patients over a given time frame (most often, 1 year):This evaluation criterion is illustrated in figure 3a. It was used in Bonhoeffer’s initial exploration of antibiotic protocols [32] (albeit with a constant offset). Under reasonable assumptions, it can be shown that maximizing X is equivalent to minimizing patient fatalities and minimizing total hospitalization time (appendix A).

Figure 3

Comparison of different optimization criteria. Here, we consider some hypothetical hospital for which the number of uninfected patients X(t) bounces around at high levels before the introduction of a small number of multiresistant infections at time . At this stage, multi-resistance increases to high prevalence (T1/2), while the number of uninfected patients drops, oscillating around a low equilibrium value, , for the rest of time. (a) the X365 optimization criterion attempts to maximize the total number of uninfected patients (area under the curve) up until a given time (usually t = 365.) (b) The T1/2 optimization criterion seeks to maximize the time taken until multiresistance takes over half the population, R = 0.5. (c), X aims to maximize the number of uninfected patients up until T1/2. (d) X maximizes the gain in uninfected patients relative to the multiresistance equilibrium, prior to multiresistance. Rather than integrating over X(t) we instead calculate the X equilibrium both before () and after () the introduction of multiresistance. While this criterion may seem unusual, one step removed from a direct integration of X, unlike the other criteria it does not depend on initial conditions.

As an evaluation criterion, X365 runs into two difficulties: firstly, it is explicitly 'blind’ to the time of multiresistance introduction. Secondly, selection of different time windows or initial conditions may lead to different results; long time horizons emphasize the eventual steady state, while shorter time horizons are informed by initial transient behaviour. It is not always clear which time window is most appropriate.

If our interest is primarily in the arrival dynamics of multi-resistant bacteria, then we may instead attempt to maximize T1/2, the time at which half of all infections are doubly resistant (R). While interesting and meaningful from an evolutionary standpoint, T1/2 makes a poor optimization criterion; maximization of T1/2 generically leads to withholding all antibiotic use, thus delaying the takeover of doubly resistant mutants indefinitely. This optimization criterion is illustrated in figure 3b.

One possible balance between these two conflicting goals is to instead maximizeThis is illustrated in figure 3c. While providing some balance between maximizing health and time, it is doubtful that X provides the correct balance between X365 and T1/2. In effect, X assumes that we will have no healthy patients from T1/2 onwards. It seems unlikely this was the intent of authors using this criterion, but it is the implicit result.

Much like T1/2, X is liable to recommend non-treatment, as integrating even small numbers over an infinite time window gives ‘optimal’ results. Also, as mentioned by Uecker & Bonhoeffer [41], the criterion completely ignores the effects of a particular epidemic protocol on the time after the emergence of R. Because X depends on the exact time course of X(t), there is also the risk that it will give conflicting answers for systems with different initial conditions. Because the exact initial conditions for a particular hospital are not knowable in advance, it would be preferable to avoid such sensitivity.

All three of these difficulties can be avoided by instead considering the novel optimization criteria (as illustrated in figure 3d)Here, indicates the long-term average of X prior to multiresistance, is the long-term average after multiresistance occurs and denotes the first time that the multiresistant population reaches some low level ε, potentially set so that ε represents a single multiresistant infection. The use of long-term averages removes the effects of initial conditions and makes X terms more analytically accessible. In the case of cycling, these averages are taken over the course of one entire cycle; in the case of ‘static’ protocols, these averages are the equilibrium values of X.

X can be thought of as being ‘formally equivalent’ to optimization over the integral , albeit with the ‘constant’ subtracted off so as to render the optimization criteria finite. In some sense, our goal is not to maximize the total number of healthy patients over some time window (which can be manipulated via manipulation of the time window), but instead the increase in health due to the use of antibiotics, over all time.

Subtracting off also serves to forbid certain degenerate strategies where . These strategies represent protocols which give worse patient outcomes than the threat of multiresistance itself (for example, never using antibiotics), but delay multiresistance indefinitely (). Protocols of this kind tend to give ‘infinite value’ according to time focused criteria such as X and T1/2. By contrast, X assigns a negative score to such protocols (). We will discuss this detail more thoroughly in §4. See figures 7 and 8.

Figure 7

(a,b) The expected values of X365 for high and low mutation/importation rate ν. Each line indicates a different channel for introduction of multiresistance. (c,d) Expected values of T1/2 and X. As might be expected, both optimality criteria tend to infinity as χ nears 1 or 0, indicating that for sufficiently extreme values, R will never account for half the infected population. This applies regardless of the multiresistance introduction channel. (e) X for a variety of χ values. (f) X rescaled such that max (X) = 1. This rescaling makes maxima more clearly identifiable, and is as mathematically valid as any other scaling, seeing as comparing ν between different multiresistance introduction methods is inherently meaningless in the context of the current model.

Figure 8

The same as figure 7, except here we optimize for different cycle times as opposed to different mixing ratios. (c,d) Once again we see that T1/2 and X are maximized for extreme values (in this case, infinity cycle time). In (f), we observe that regardless of the multiresistance arrival channel, the saturation time tsat gives close to optimal results.

It is important to note that no assumption is made that and are achieved using the same management protocols, hence (for example) it is possible that a cycling protocol may be used prior to multiresistance, and a mixing approach afterward. X optimizes on the assumption that we pick the best possible protocol post multiresistance introduction.

is chosen over T1/2 as our stopping criterion for technical reasons. See appendix B.4 for details, as well as a comparison between X and the G1/2 optimization criterion put forward by Bonhoeffer [32].

See figure 3 for a schematic illustration of all optimization criteria. Optimization criteria are also summarized in table 1.

Mean X values

In order to make sense of X (and other optimization criteria), it will prove helpful to calculate both and . Rather than calculate numeric integrals of X for a variety of initial conditions and parameter values (as has been done in previous papers [33,35,36]), our goal in what follows is to find asymptotic approximations for a variety of parameter regimes, making use of different simplifying assumptions in each case. See figure 4 for sample trajectories in each of these regimes.

Figure 4

General behaviour of the system when using mixing (a), combination therapy (b), or cycling (c,d). For the two ‘static’ strategies, (a,b), random initial conditions rapidly approach equilibrium for S and X. In the mixing case (a), R and R approach equilibrium slowly; for the parameters considered here, they become equal in limit as t → ∞; differences are purely based on initial conditions. For cycling (c,d), background colour indicates the antibiotic currently in use. We consider two distinct regimes: ‘fast cycling’ (c), in which equilibrium is never reached, and ‘slow cycling’ (d), in which R, R approach equilibrium values and are stabilized by importation from the community before drug switching occurs.

Let us start by considering the case of combination therapy prior to the introduction of multi-resistant bacteria. In this case, assuming antibiotics are effective at keeping down infection (τ + μ) ≫ β X for each β, all infected compartments can be well approximated by the balance between immigration and recoveryThe total population of the ward is given by , and hence the long-term equilibrium of X is well approximated by

In the case of mixing (prior to multi-resistance), we use the 7-box Uecker model (figure 2). In order to calculate the long-term equilibrium behaviour, we first calculate what fraction of cases are receiving effective treatment for each strain. Effective treatment ratios can be shown to be equal toand

This in turn leads to three different approximations of X at equilibrium, depending on which resistant strain is dominantThe true equilibrium value for is well approximated by the smallest of these three values. If is smallest, this indicates that R is the dominant infection strain, the current limiting factor on improved health. If gives the smallest value (usually for high treatment correction values, q), this indicates that infection is dominated by disease importation; treatment within the hospital is close to optimal and ABR within the hospital is dominated by importation from the community. This situation might arise for example in situations where we have high levels of testing for antibiotics (such as the Dutch ‘search and destroy’ ABR policy [49]), or in situations where community levels of ABR are very high.

Illustrations of these results are given in figure 5. Derivation of these results can be found in appendix B.2.

Figure 5

Equilibrium values of uninfected individuals when using a mixing protocol. (a) As q (the drug correction rate) is increased, increases approximately linearly (following the equilibrium), before approaching a maximal value at the equilibrium. (b) As drug selection probability χ is varied, the equilibrium population increases to a maximum and then decreases, passing from a limited regime to a limited regime. Here, we select drastically different β, β, so as to illustrate the effects of asymmetric infection rates. In practice, these rates can be expected to be approximately equal, and χ = 0.5 gives close to optimal results. For both figures, the solid blue line represents numeric estimation of X at equilibrium, calculated using Newton’s method. The black dashed lines are the upper bounds, as given by , , above. The black dotted line (just below the dashed line, left panel) indicates a somewhat tighter upper bound for ; this improved approximation gives only modest improvements, and requires significantly more algebra, see appendix B.2 for details.

The next case to consider is the cycling case. We are interested in the long-term mean value of X averaged over precisely one cycle length: , where T is the length of a single cycle (we consider t0 → ∞ so as to ignore transient effects that might be caused by initial conditions. Alternatively, we might solve to find X(t) on the interval [0, T] with periodic boundary conditions). For the time being, we consider the symmetric case, in which A and B have identical migration and treatment parameters (m = m, β = β) and equal cycle time. Detailed calculations applicable to both the symmetric and asymmetric case are provided in appendix B.3.

For cycling, there exist two major parameter regimes: a ‘fast cycling’ regime in which the equilibrium is never reached and a ‘slow cycling’ regime in which the system approaches its long-term equilibrium with each cycle (figure 4c,d). For sufficiently fast cycling, it can be shown thateach antibiotic is used half the time, and hence, each applies at 50% effectiveness. For slow cycling, the long-term average approachesHere, C represents the transient ‘spike’ in X immediately following the change in treatment, while resistance to the new drug is rare (figure 4). The effects of the spike are diluted across the length of a cycle. represents the equilibrium value of X approached over the course of an arbitrarily long cycle, once resistance to the new drug has taken hold. Each drug has less than 50% effectiveness, because slow cycling results in resistance becoming ubiquitous in the population with each cycle; each drug spends most of its time treating the bacterial strain it is least effective against. So long as resistance importation is rare (m, m ≪ R, R) both and C can be calculatedAs previously, σ is the total population. and are the equilibrium population size for susceptible and resistant strains, assuming appropriate treatment (so, R being treated with B, or R being treated with A).

For all T, is best approximated by the minimum of , , see figure 6. This approximation is generally fairly tight, except in the boundary region where . We refer to the boundary between fast and slow cycling as the ‘saturation time’, denoted tsat. Saturation time can be calculated by setting , and solving

Figure 6

Comparison of as calculated numerically, versus the approximation . Analytic results provide an accurate approximation of for most T values, except in a small window on either side of the ‘saturation’ time tsat.

The final parameter regime to consider is the equilibrium in the presence of multiresistant bacteria, that is to say R > 0. In this case, it is easy to show that R must increase until

This result applies regardless of which treatment protocol is employed. Strategies that would receive higher in the absence of multiresistant infections will instead approach an equilibrium or cycle with mean value close to . Protocols with worst health outcomes than the best possible response to multiresistance (), such as withholding antibiotics entirely, will retain their low X value, and will result in the R population decaying at a rate proportional to (see appendix B.4 for details).

Introduction of multiresistance

Analytic calculations of before and after the introduction of the doubly resistant strain provide some measure of the relative ranking of different ABR protocols. In order to evaluate however, we must also consider when double resistant bacteria are introduced, that is to say, we must estimate . Because X is linear in , it is sufficient to estimate the expected value of .

Doubly resistant strains are introduced to a hospital via a number of different channels, each with its own unique rate constant. For example, if a doubly resistant strain is imported from outside the hospital at some constant rate Mimport = m, then the expected introduction time of doubly resistant infection is constant, and independent of the antibiotic management protocol currently in use in the focal hospital (though potentially dependent on other hospital or inter-hospital policies [55]).

Multiresistant strains can also arise via de novo mutations. Such mutations can be modelled as either dependent on selective pressure () [34] or entirely independent of selective pressure () (one of multiple possibilities originally considered by Bonhoeffer et al. [32]). Finally, multiresistance may also occur via horizontal gene transfer between existing strains; the case was considered by Bergstrom et al. [33], where the assumed multiresistance is produced at a rate proportional to MHGT = νHGT(R × R). More detailed models of HGT, based on within host dynamics, have been studied (for example Roberts et al.'s recent work [56]). For the time being, and for the sake of simplicity, we follow in the footsteps of Bergstrom et al. The relative importance of each of these four channels can vary from pathogen to pathogen, and will also depend critically on antibiotic stewardship policy across the surrounding region; careful stewardship may mean that Mimport is minimized and we are mainly interested in de novo mutations or gene transfer. In a region with high ABR prevalence Mimport may dominate.

While each of the above are reasonable assumptions they are all, in some sense, ‘cartoon’ approximations of exceptionally complex processes. Future investigation into the various sources of multiresistant infection may well suggest improvements upon the above terms, or even introduce new terms representing previously ignored channels such as horizontal gene transfer from a patient’s commensule bacteria [57]. At the very least, it would be useful to determine the relative contribution of each channel in clinical practice. Such questions, however, are far beyond the remit of this simple mathematical analysis. For the time being, we treat each of the above M as given and assume that in any given circumstance one channel of multiresistance dominates all others.

Our goal in what follows is not to make any specific or universally applicable policy recommendations, but instead to examine the types of recommendation made by various optimization criteria under the influence of different M. We are in some sense evaluating the optimization criteria themselves, rather than the management protocols on which they act. Rather than attempt a full of exploration of the entire parameter space, we focus on two straightforward test cases in order to illustrate the types of behaviour generally observed.

Comparison of optimization criteria: mixing

In what follows, we make use of Uecker’s 7-box model, with parameter values β = 1, β = β = 0.99, β = 0.98, q = 0. We are interested in the recommendations made by the four optimization criteria X365, T1/2, X and X with respect to the mixing rate parameter χ. For the sake of breaking symmetry slightly we set m = 2 m. Results in the q > 0 case are explored in appendix D. A full description of all parameter values and simulation methods is made available via github [58], and is saved in the file ‘EvaluatingAllOptimalMetricsMixing.m’.

In all cases, we assume that R(0) = 0 and that multiresistant mutants appear according to a Poisson process with rates proportional to Mimport, Mbase, Mselect or MHGT. Results for each of the optimization criteria are presented in figure 7. Both X365 and X have local maxima near χ = 1/2 for Mimport, Mbase, Mselect. For MHGT, these two optimization criteria have sharp local minima near χ = 1/2, with their maxima off centred (to differing degrees). This makes sense: if resistance is primarily formed via horizontal gene transfer, than the primary goal of any management protocol is to avoid regions of treatment space where R and R coexist. X365 and X differ in how far from χ = 1/2 one must move in order to reach optimal results.

By contrast, X and T1/2 both recommend extreme values of χ, tending to infinity precisely in those areas where X < 0. Once again, this matches expectation: if then R cannot increase in prevalence, and hence T1/2 = ∞.

When dealing with the optimality criteria X, X and T1/2, ν acts primarily as a scaling constant (with a constant offset in the case of X and T1/2, caused by the delay between T1/2 and ). In the case of X365, by contrast, the value of ν acts to determine the competition between multiresistance arriving and the year ending, leading to qualitatively different results depending on the exact value of ν (figure 7a,b). Fortunately, these changes have at most modest effects on the optimal mixing ratio, even for X365. Hence, precise knowledge of ν is not crucial.

The important lessons from figure 7 in terms of criteria comparison are that X and T1/2 are virtually indistinguishable, and can be seen to optimize in almost precisely the opposite direction to X365 and X. X365 and X are broadly similar in their overall behaviours, though not always in their precise recommendations.

Comparison of optimization criteria: cycling

We next consider the implications of different optimization criteria in the context of antibiotic cycling. In this case, the 5-box model is used, with parameter values β = 1, β = β = 0.99, β = 0.98. Rather than attempt to explore the entire parameter space of possible cycling protocols, we here assume symmetric cycling times, and equal importation rates (T = T, m = m), see appendix D for comparisons in the asymmetric case. We are interested in determining the expected value of X365, T1/2, X and X for a variety of cycling times T. Expectation is taken over all possible introduction times of the multiresistant mutant, and also over the ‘initial phase’ (how far through the cycle the system is at t = 0). Initial phase is selected uniformly at random between 0 and 2T. Full code is available via github [58], and is saved in the file ‘EvaluatingAllOptimal_cycling5box.m’. See figure 8 for results.

Once again, X and X give conflicting recommendations. Much like the mixing case, X and T1/2 both recommend extreme parameter values: both metrics increase monotonically with cycle time T and tend to infinity as T → 3957.86; this is the smallest value for which . By contrast, both X365 and X are maximized for intermediate values of T. For X, optimal cycle times for all multiresistant introduction channels cluster around tsat. This makes sense—tsat denotes the largest T value that can be used without suffering reductions in . Any M that is minimized by increasing cycle time T can be increased up to tsat at no cost. Any increase beyond tsat inevitably comes at the cost of a reduction in . Because the transition between and at tsat is not sharp, the exact location of the maximum of X varies slightly depending on the source of multiresistance. For Mbase and Mimport, X is maximized just below tsat, for MHGT and Mselect, the cost of switching antibiotic is higher, and X is maximized for T slightly larger than tsat. This clustering of optimal results around tsat is stable to variations in parameter values; see appendix C.

Qualitatively, X365 gives results that are similar to X in some ways, but not identical. Overall, X365 can be described as ‘lumpier’; this higher complexity results from the interaction between three different timescales: the timescale of cycling, the timescale of ABR introduction, and the timescale of the integral (one year). When ν values are very small, the probability of a mutation occurring within 365 days becomes small. In this case, X365 approaches . By contrast, when ν values are large, mutation occurs within one or two cycles. In this case, the approximation is no longer appropriate and the exact phase of cycling at the start of the simulation has a significant impact: for example, if using MHGT then multiresistance most commonly arises during the time of drug switching. Whether t = 0 occurs before or after a switch can have significant impact on the time until multiresistance arrival, and hence on the overall shape of X365. The large ν regime for X365 draws attention to transient effects caused by initial conditions. These time-scale effects are ignored by X, which considers only long-term averages of M, and has no predefined ‘end point’ at 1 year.

Comparison of optimization criteria: combination therapy

The final antibiotic management protocol that we consider is combination therapy, in which both antibiotic treatment options are used simultaneously for all infections. Unlike the previous two cases, where we have a parameter value to vary over, combination therapy has no such parameter values: it is, in essence a single protocol, as opposed to a wide family of protocols. As such, there is no good way of comparing combination therapy to itself, and instead the protocol must be compared to the best results from the previous protocols. The results below are based on solving equation (3.4) and comparing to the optimal results for mixing and cycling. Full code for this work can be found on github [58] included as a special case in the file ‘EvaluatingAllOptimalMetricsMixing.m’. The parameter values assumed were μ = 1/5, β = 1, β = β = 0.99, β = 0.98, γ = 1/10, γ + τ = 1/2.5, other import/export parameters can be found in the file itself, as needed.

Using either cycling or mixing, X and T1/2 can be pushed towards infinity for suitably extreme parameter values (χ → 0 or 1 in the case of mixing, T → ∞ in the case of cycling). Combination therapy allows no such ‘infinite optimization’ for either of these protocols, and hence, for any protocol focused on the dominance time of multiresistance, must be considered strictly worse. This is in line with results by Obolski & Hadany [34], who rank protocols based on the emergence time of multiresistance, and conclude that cycling is preferable to combination therapy.

For X365, with either large or small ν values, combination therapy is of the order of 2-3 times better than optimal mixing and optimal cycling for Mbase, Mselect and MHGT (cycling beats mixing for Mbase, Mimport, but is inferior for MHGT, though in most cases, differences are marginal). For Mimport, we find that optimal cycling is superior to combination therapy, which is superior to optimal mixing, regardless of ν.

For X however, we find that combination therapy dominates all other strategies, regardless of M. For Mimport, combination therapy gives X three times larger than both optimal cycling and optimal mixing. For Mselect, we find a approximately 190-fold increase, for Mbase, combination therapy is approximately 500 times better than optimal (symmetric) cycling, and approximately 700 times better than optimal mixing. In the case of MHGT, combination therapy gives X values more than 1000 times greater than optimal mixing, which is in turn more than 10 times greater than optimal symmetric cycling (this difference is reduced for asymmetric cycling, but in this case, the optimal cycling times tend to zero, indicating that mixing is the superior strategy, see appendix D for details). The large ratios involved here should be read with a degree of caution: we do not claim that combination therapy is hundreds (or thousands) of times ‘better’ than mixing or cycling strategies. By its nature X takes values closer to 0 than many of the other evaluation criteria, hence exaggerating the ratio between different values. X365 may be preferable for more physically meaningful comparisons between combination therapy and other protocols. That said, for models considered here it is clear that X unequivocally favours combination therapy over all other protocols.

Discussion

Antibiotic resistance, and the proliferation of multi-resistant bacteria, pose significant challenges to modern healthcare systems, threatening to roll back the past century of antibiotic research [16]. Antibiotic management protocols are designed with the goal of improving patient outcomes while preventing (as much as possible) increases in resistance. The question of how best to represent ‘good outcomes’ mathematically runs into certain difficulties: individual optimization criteria often run at cross purposes and in many cases entirely contrary to one another.

Based on our explorations in §4, it seems likely that criteria intended primarily to delay the arrival of multi-resistance (T1/2 and X) should be avoided in most circumstances. These criteria tend to infinity precisely in those regions where health outcomes are worse than the long-term impact of multiresistance itself. This may be appropriate in certain cases: when dealing with mild, short-term illness, antibiotic stewardship may be prioritized over immediate recovery [59]. This is not generically the case considered in this article however, where our focus has been antibiotic policy for preliminary antibiotic allocation for inpatients at a hospital (prior to more detailed ABR testing). With this in mind, it would appear that time maximizing optimality criteria are most often actively harmful to the patient population, both in the short and long term. It is precisely these time maximizing optimality criteria that recommend against combination therapy [36]. In all other cases, when combination therapy is considered, it is found to be superior to both cycling and mixing based protocols. Whether or not this rather strong result holds up when considering additional factors such as monetary cost or bystander selection is currently unknown.

In order to balance the value of delayed multiresistance with improved health outcomes, we construct the novel optimization criterion, ; this is in some sense equivalent to optimizing , albeit with the ‘constant’ subtracted off so as to avoid infinities. By definition, the criterion only gives positive values for protocols that improve patient outcomes relative to a ward dominated by multiresistant bacteria.

Asymptotic arguments allow us to calculate the mean number of uninfected patients, , for a variety of cases, both pre- and post-multiresistance arrival (§3). can be shown to be independent of management protocol, while is protocol dependent. In almost all approximations of , the infection rate β is a key parameter. While we have generally discussed differences in β as being metabolic costs of resistance, it is important to note that these costs can be imposed ‘artificially’ through targeted isolation of infected individuals. This is suggestive of the critical importance of such clinical measures as improved hygiene practices and rapid diagnostics and isolation of ABR cases [48,49,60].

In all cases, we find that our results are sensitive to ABR importation rate m, m and m, that is to say, the prevalence of ABR in the community. When using cycling, we find that optimal cycle times for scale with tsat, the so-called ‘saturation time’ above which which rapidly decays (figure 8). The exact position of optimal cycling relative to tsat depends on the source of multiresistant infection (horizontal gene transfer, spontaneous de novo mutation, or selective de novo mutation). Because tsat depends on the prevalence of ABR in the community, knowledge of the local community may be crucial for selecting optimal cycle times.

Recommendations for optimal mixing depend on the means of multiresistance introduction: balanced 50:50 mixing gives good results when multiresistance is either imported, or generated through de novo mutation and poor results if multiresistance arises via horizontal gene transfer (figure 7). Further empirical work will be needed in order to determine which of these channels is most significant.

While the discovery of tsat, and its use in estimating optimal cycling times is a nice result, there are (inevitably) a number of caveats, conditions and stones still left unturned.

Firstly, it is worth reiterating that we have deliberately based our work here on a classical model [32]. This was done to allow cleaner comparison, and give readers a baseline that we hope they are familiar with. With that said, science moves on. Our modelling tools and understanding of biology have improved since Bonhoeffer’s publication. The significant impact of phenomena such as bystander selection [61] are better understood, and there has been a more recent push toward integrating both within host and between host models [62]. Similarly, in the area of clinical policy, diagnostic-informed strategies such as the ‘search and destroy’ ABR policy implemented in the Netherlands [48,49,60] have been mentioned only briefly in this article. Such diagnostic informed policies can be hugely powerful, and were neglected here simply because it was our intention to model pre-diagnostic ABR policy (in line with the literature we wished to comment upon).

As a more mathematical caveat, many of the asymptotic results here are made on the assumption that both m and m are rather small; ABR spread is dominated by infection within the hospital (nonsocomal infection). Outside of this parameter range, the results presented here may be less relevant. The second mathematical limitation in the present research is the assumption of continuously varying population; given the relatively small size of a hospital (dozens to hundreds of individuals), this continuity assumption is at best suspicious, especially when many key dynamics of the system occurring when R(t), R(t) ≪ 1. It would be interesting to explore these results in the stochastic context. It seems likely that some analogue to both and T can be determined, though it is far from clear that will be independent of ABR protocol in the stochastic case. Similarly, the addition of a third antibiotic is also predicted to change the behaviour of the system, and the appropriate definition of . Despite these complications, we would still posit that some optimization criteria conceptually equivalent to X are likely to prove useful in all of these cases.

While was found analytically for all management protocols, calculation of mutant arrival rates relied on numeric solutions of ODEs, and no analytic approximation looks forthcoming. This problem is further exacerbated by the fact that it is not even clear which M function is appropriate, and the possibility of horizontal gene transfer between native commensule bacteria carrying resistance genes and invasive pathogenic bacteria [63-65] raises the possibility that none of the M functions explored here give reliable results. Determination of which M function should guide selection of cycle time is an empirical question rather than a mathematical one. We do note, however, that while seldom optimal, T = tsat is considered ‘good’ for all channels of multiresistance considered here.

Throughout the literature, numerous antibiotic deployment protocols have been proposed, each with the conflicting goals of maximizing patient health and maintaining antibiotic effectiveness. Our focus throughout this article has been to compare these protocols and (more crucially) to compare the various optimization criteria used to rank them. We find criteria that are overly focused on antibiotic stewardship (T1/2, X) tend to recommend patient outcomes which are worse than the long-term outcome of multiresistance, and hence are harmful to patients both in the short and long term. For this reason, we recommend against the use of such optimization criteria. We find that optimization criteria which are more patient-centric (X365, X) generally recommend combination therapy as the best method of preventing the creation of multiresistance inside the hospital, but may occasionally recommend cycling if multiresistance is primarily introduced from the community.