The Dosing Schedule Decides Whether Resistance Wins

The Poster in the Waiting Room

Every one of us grew up with the same instruction. Finish the course. Take the last pill even when you feel fine, because stopping early is what breeds resistant bacteria. It is on posters in surgeries and on the leaflet folded inside the box, and most of us had said it out loud to somebody else without ever having checked it.

So we checked it, in the only way a school club with no laboratory can check anything. We wrote down a population model, gave it a drug concentration that rises at each dose and decays in between, put a sensitive strain and a resistant strain in it, and ran the thing 178,002 times.

The answer was not the one on the poster, and it was not the opposite of the one on the poster either. It was a third thing. In our model the height of the concentration peak decides almost the whole outcome, and the length of the course stops mattering as soon as the drug has finished its work. Missing doses does real damage, for a reason that has very little to do with how much drug you end up swallowing.

Two warnings before anything else, and we mean both of them. First, this is a simulation. No bacterium, no patient and no clinical record was involved at any point. When we write that a schedule "cured" 86.7% of cases, we mean that 1734 of 2000 simulated populations in a computer went to zero. Second, and more important, nothing here is advice about medicine. A person deciding when to stop taking an antibiotic is facing a problem that includes their immune system, the site of their infection, the specific organism and a doctor who has seen the patient. Our model contains none of those. It is a study of a mathematical object that resembles an infection in a few respects and differs from one in many. We will keep saying so.

What Is Actually in the Model

The population contains two strains of one species. The sensitive strain S grows at \(\psi^{\max}_S = 1.0\) per hour when there is no drug, which is a doubling time of about 42 minutes. The resistant strain R carries a fitness cost \(c = 0.10\), so its drug-free maximum is \(\psi^{\max}_R = 0.90\) per hour. Costs in that range are typical of chromosomal resistance mutations before compensatory evolution erodes them [9].

Drug acts through the sigmoid pharmacodynamic function of Regoes and colleagues [3], which fits measured net growth rates for five antibiotic classes and is the standard way to write killing as a function of concentration [2][4]. For a strain with maximum net growth \(\psi^{\max}\), maximum kill rate \(\psi^{\min}\) and MIC \(m\):

We used \(\psi^{\min} = -4.0\) per hour for both strains and a Hill coefficient \(\kappa = 1.5\). The function is built so that \(\psi(0) = \psi^{\max}\), \(\psi(\infty) = \psi^{\min}\) and, the part that matters, \(\psi(m) = 0\) exactly. That last identity is what an MIC means: the concentration at which net growth is zero. A little algebra collapses the whole expression to

which is the form in the code. The sensitive strain has \(m_S = 1\) mg/L and the resistant strain \(m_R = 16\) mg/L, a sixteen-fold shift of the kind a single efflux or target mutation can produce.

Pharmacokinetics are deliberately the simplest thing that still has the right shape. One compartment with instantaneous input, and first-order elimination at a half-life of 3.0 hours, so \(k_{\mathrm{elim}} = \ln 2 / 3 = 0.231049\) per hour. Each dose that is actually taken adds \(D\) mg/L to the concentration and the concentration then decays. A regimen is therefore two numbers, the dose \(D\) and the interval \(\tau\), plus how many doses and how reliably they are swallowed.

Demography is a stochastic birth and death process on integer cell counts, advanced by tau-leaping on a fixed step of 0.02 h [16]. Births are Poisson draws at rate \(\psi^{\max}_i\), independent of drug. Deaths are binomial draws on the standing population at rate

so that the net growth of strain \(i\) is \(\psi_i(C) - \psi^{\max}_S N_{\mathrm{tot}}/K\). With no drug and the population at carrying capacity \(K = 10^9\), the sensitive strain sits exactly at equilibrium and the resistant strain declines at rate \(c\), which is what a fitness cost is supposed to do. Drawing deaths as binomials rather than Poissons means a population can never be pushed below zero by a leap that is slightly too large.

Resistance arises only by mutation, at \(\mu = 10^{-9}\) per sensitive birth. There is no back mutation and no horizontal transfer. Each run begins with 48 hours of drug-free burn-in from \(N_S = K\) and \(N_R = 0\), which lets standing variation build up stochastically instead of being assigned by hand, and then the first dose lands.

Nothing Gets Reported Until the Window Checks Out

A simulation that has never been checked against something outside itself is a picture, not a measurement. So the program runs seven validations before it produces a single result, and it prints our value beside the accepted value with the difference, every time.

The first two concern the mutant selection window, the concentration band in which resistant mutants are selectively amplified [5][6][7]. Drlica and Zhao define its upper edge as the mutant prevention concentration and its lower edge, in the classical statement, as the MIC of the susceptible strain. In our model both edges are exactly where a strain's net growth crosses zero, so bisecting on \(\psi_S\) and \(\psi_R\) should return the two MIC values we typed in. It does, to 2.9 parts in \(10^{15}\).

The wider window is more interesting. Selection for resistance does not wait for the sensitive strain to stop growing; it begins as soon as \(\psi_R\) exceeds \(\psi_S\), which happens below \(m_S\) whenever the fitness cost is small. Setting the two expressions equal and substituting \(X = C^{\kappa}\) turns the equation linear, and it solves:

with \(p = \psi^{\max}_S\), \(q = \psi^{\max}_R\) and \(C_{\mathrm{sel}} = X_{\mathrm{sel}}^{1/\kappa}\). One of us derived that on paper before the code was written, which is the only order in which such a check is worth anything. The bisection and the closed form agree to 2.7 parts in \(10^{15}\).

Two of those deserve a sentence each. The zero-cost check asks whether the model can make the resistant strain decline for a reason it should not have. Run deterministically with \(c = 0\) and no drug for 500 hours, the smallest change in \(N_R\) over any step was exactly zero and the smallest net growth rate was exactly zero, so it never declines. Put the real cost back and the same run takes \(N_R\) from \(1.0 \times 10^3\) to \(1.9 \times 10^{-18}\), a decay rate over the last hundred hours of 0.100000 per hour against the 0.100000 we typed in.

The zero-mutation check is the one that surprised people. Take a schedule that is hopeless in the full model, 8 mg/L every 24 hours for a week, which cures only 8.90% of replicates. Set \(\mu = 0\) and start with no resistant cells. It now cures 2000 out of 2000. The schedule is not failing because it is too weak to kill bacteria. It is failing because it is strong enough to clear the sensitive strain and too weak to touch the resistant one.

Notes From the Club Table

What One Missed Dose Buys, in Arithmetic

This section has no simulation in it. It is the pharmacokinetics and the pharmacodynamic function integrated directly, and it explains almost everything that follows.

Start with the concentration. With doses of \(D\) every \(\tau\) hours and first-order elimination, the steady-state peak is the geometric series \(C_{\max} = D/(1 - e^{-k\tau})\). For \(D = 32\) mg/L every 8 h with \(k = 0.231049\):

So on this schedule the drug spends 47% of every interval above the concentration at which the resistant strain cannot grow, and it never falls anywhere near the bottom of the selective window before the next dose arrives. Now integrate the resistant strain's net growth over one interval, \(\exp \int_0^{\tau} \psi_R(C(t))\,dt\), and over one interval that has been doubled by a missed dose:

Read the two blocks against each other. On the high schedule a resistant cell ends every interval slightly worse off than it started, because 0.5598 is less than one. It is losing ground on the drug's own terms. One missed dose turns that 0.5598 into 386, a swing of a factor of 690, and roughly 1.93 such gaps are enough to take the 9.9 resistant cells the model starts with up to a million. That is where the adherence result in §9 comes from.

On the low schedule the resistant strain multiplies by 475 every single interval whether or not anybody misses anything. Perfect adherence on a schedule like that is perfect adherence to a selection experiment.

Twenty-Five Schedules, Run at Perfect Adherence

The main grid crosses five dose sizes with five dosing intervals, gives every cell a seven-day course and perfect adherence, and runs 2000 replicates in each. That is 50,000 simulated courses. We report the standard error from the replicate counts and also a 95% Wilson score interval, because several cells come back at exactly zero or exactly one, where the plain normal interval has zero width and says something false.

Read down the fourth column and the usual intuition fails. The cell at 8 mg/L every 4 hours pours 1454.2 mg/L·h of drug into the model and cures 11.05% of replicates. The cell at 32 mg/L every 24 hours delivers 969.5, a third less, and cures 68.70%. More drug, worse outcome, by a wide margin in both directions.

The column that does order the table is fT>MIC_R. Clinical pharmacology has spent decades asking which summary of a concentration curve predicts outcome, and the three standard candidates are the fraction of time above MIC, the peak relative to MIC, and the daily area under the curve relative to MIC [4]. We computed all of them for our 25 regimens and ranked them against the measured cure rate:

Peak height wins. Because the peak divided by MIC_S and the peak divided by MIC_R differ by a constant factor of 16, a rank correlation cannot tell them apart and both come out at exactly +0.9809; the ranking says peak height matters and is silent about which MIC to divide by. Time above the sensitive MIC, the index that gets quoted most often for beta-lactams, is nearly useless here at +0.6078, because almost every schedule in the grid clears that bar. And time spent inside the selection window correlates negatively with cure, which is the mutant selection window hypothesis written as a number [6][7].

There is a second, blunter way to see what is going on. Split all 50,000 replicates by whether a single resistant cell existed at the moment the first dose was taken. The 3479 replicates that started clean were cured 67.98% of the time. The 46,521 that started with at least one resistant cell were cured 25.12% of the time. The odds ratio is 6.33. Almost everything the treatment schedule does, it does to a resistant population that was already there before the first pill.

Same Drug, Different Shape

If the peak is what matters, then two schedules delivering identical total drug should give different answers as long as they concentrate that drug differently. We built the ladder to test it. Dosing interval fixed at 8 hours, the product of dose size and dose count held at 288 mg/L, so every arm delivers an AUC of exactly 1246.5 mg/L·h. Only the shape changes.

The gap between the top row and the bottom row is 0.9110 in probability of cure, which is 143 standard errors, on the same amount of drug. Spread a fixed supply of antibiotic thin enough and you convert a cure into a selection experiment. That is the aggressive versus moderate question that the evolutionary medicine literature has been arguing about for fifteen years [11][12][13], and on this model, at this inoculum, aggressive wins by a distance.

But the matched ladder changes two things at once. Going from row one to row seven the course gets longer and the dose gets smaller, and we cannot tell from those numbers alone which one did the work. So we ran a second ladder holding the dose fixed at 32 mg/L every 8 hours and varying only the number of days.

Hold the dose and the picture inverts. One day of 32 mg/L every 8 hours cures 86.70% ± 0.76%. Two days cures 96.70% ± 0.40%. Three days cures 99.55% ± 0.15%, five days 99.90% ± 0.07%, and seven, ten and fourteen days all cure every one of 2000 replicates. So course length does help. It helps steeply for the first two days and then stops, and the fourteen-day arm buys its last 0.45 percentage point of cure with 4570 extra mg/L·h of drug.

Six Hours

Here is the number the club kept coming back to.

In the fixed-dose ladder, the median time to clearance was 5.88 h in the one-day arm, 6.02 h in the two-day arm, 6.04 h at three days, 6.06 h at five, 6.00 h at seven, 6.04 h at ten and 6.06 h at fourteen. Seven different course lengths, spanning a fourteen-fold range of total drug, and the infection ends at the same moment in all of them.

It has to. The drug does not know how many doses are left in the packet. Once the sensitive population has gone from \(10^9\) to zero, which at a kill rate near 3.7 per hour takes about six hours, every remaining dose is being taken by a patient who no longer has that infection. In this model those doses are not neutral. They are the period during which any surviving resistant cell has the field to itself and a drug concentration that is often sitting in its selection window.

That is the honest version of the disagreement with the poster. The poster is right that stopping too early is dangerous, and our one-day arm shows it: 13.30% of replicates ended in resistance because treatment stopped while sensitive cells were still there. The poster is silent about the other end, where continuing past clearance is also dangerous, and our matched ladder shows that one too. The safe zone in this model is narrow and it is early. Clinical trials of shorter courses have been finding something with a similar shape for a decade, without needing our model to tell them [17][18].

A Gap Is Not the Same as a Smaller Dose

Now the missed doses. The reference is 32 mg/L every 8 hours, the schedule that cured 2000 out of 2000 when taken properly, and each scheduled dose is now taken independently with probability \(a\).

Missing one dose in five costs 14.05 percentage points of cure, 18.1 standard errors. That is a real effect and it goes in the direction the poster predicts. The three-day course behaves almost identically, 0.9945 down to 0.8555 at the same adherence, so the damage is done by the gaps and not by the shortening of an already short course.

Then we ran the control, and this is where one of us said out loud that the model must be broken. Give every dose, but shrink each one by the factor \(a\), so that the average drug delivered matches the missed-dose arm exactly and the only difference is that there are no gaps in it.

The ragged schedule is better than the smooth one at every level, and at 50% adherence it is nearly four times better. The model is not broken. Go back to the arithmetic in §5. On a 32 mg/L schedule the first two or three doses are usually taken in full, the sensitive population is gone in six hours, and the resistant cells are being killed during the 47% of each interval when the concentration is above 16 mg/L. Cutting every dose to 25.6 mg/L lowers the peak to 30.4 mg/L and takes that fraction down; cutting to 16 mg/L leaves only 9.27% of each interval above MIC_R, and the grid row for 16 mg/L every 8 h cures 19.45%.

We want to be careful about what this does and does not say. It does not say that missing doses is harmless, because the ragged arm loses a quarter of its cures by 50% adherence. It says that in this model, the harm from a missed dose comes from the gap it opens rather than from the drug it removes, and that a schedule which never reaches a high enough peak is worse than one which reaches it intermittently. The real-world version of the control would be a person taking half a tablet every time, and nobody does that, so the comparison is a diagnostic rather than a scenario.

The Strongest Objection We Can Make

The objection is this. Our model has no immune system, so the drug does all the killing and finishes in six hours, which makes short courses look good by construction. Ankomah and Levin built a closely related within-host model that does include an innate immune response, and their conclusion was the opposite of ours: high dose and full term [14]. If the difference between their answer and ours is one missing term, our answer is an artefact and theirs is the one to believe.

We take it seriously enough that we put it in the sensitivity run. Adding a constant immune clearance rate to both strains moves the long mild arm from 0.1100 to 0.1250 at 0.10 per hour, 0.2730 at 0.30 per hour and 0.4470 at 0.50 per hour, while the short sharp arm stays at 1.0000 throughout. So immune killing helps the weak schedule considerably and never hurts the strong one, and the short arm still wins by 0.553 even when a constant immune term is killing at 0.50 per hour on top of the drug.

That is not a full answer and we should say why. A constant extra death rate is a poor model of immunity. A real innate response grows with the pathogen load, saturates, and arrives on a delay, and Ankomah and Levin's result depends on exactly those dynamics: the drug suppresses the population long enough for the immune response to build, and stopping early removes the suppression before the build is finished. Our constant term cannot represent that at all. The most we can say is that the short arm's advantage survives the crudest version of the criticism, and that testing the real version would need a different model than the one we wrote.

A second objection, from the other direction. Alexander and MacLean showed experimentally and theoretically that stochastic loss keeps most resistant lineages from ever establishing, even at concentrations inside the selection window [15]. Our model agrees with that mechanism and our numbers show it working: the 3479 replicates that began with zero standing mutants were cured 67.98% of the time on schedules where the overall rate was 25.12%. Where we part company is scale. Their argument is about single cells; our populations begin with a mean of 10.12 resistant cells already present, because \(\mu K / c\) with \(K = 10^9\) is ten. Establishment failure from a single cell is likely; establishment failure from ten independent cells over a week of favourable selection is not. Whether a real infection carries ten pre-existing mutants or none is a question about inoculum size, and it is the single assumption our conclusions are most sensitive to. See §11.

A third objection is about the numbers themselves. A Monte Carlo estimate is only as good as its sample, and a 2000-replicate cell has a standard error near 0.01 on a probability near 0.5. Panel a of Figure 5 is the audit. Taking one cell out to 20,000 replicates, the estimate settles at 0.10395 ± 0.00216 and every intermediate running value sits within 1.67 standard errors of it. The relative error shrinks as 1/√n to within 2%, which is what it should do if the replicates are independent and nothing in the sampling is correlated. That does not make the model right. It makes the numbers we report about the model reliable to the precision we quote them.

Where a Different Choice Would Have Changed the Answer

We ran the two extreme matched arms through fourteen variants of the model, 1000 replicates each, and recorded the gap in cure probability between them. The short arm won in all fourteen. The size of the win is another matter.

Three of those rows would have given a different article.

Inoculum. Drop the bacterial load from \(10^9\) to \(10^7\) cells and the fourteen-day arm cures 97.80% instead of 11.00%, and the whole finding nearly disappears: the gap falls from +0.880 to +0.022. The reason is the arithmetic of standing variation. At \(K = 10^9\) the model starts with ten resistant cells on average, at \(10^7\) it starts with 0.1, and a population that usually contains no resistant cells at all cannot have resistance selected out of it. Everything in this article is a statement about large infections. For a small one the schedule may barely matter, which is close to the point Alexander and MacLean make [15].

Size of the resistance step. We gave resistance a sixteen-fold MIC shift, which a doctor's peak concentration can still exceed. Make it thirty-two-fold and the short sharp arm collapses from 0.9900 to 0.1950, because 37.98 mg/L is no longer above MIC_R = 32 for any usable fraction of the interval. The entire high-dose strategy in this article rests on being able to reach above the resistant strain's MIC, and real dosing is bounded by toxicity, which our model has no representation of whatsoever [13]. Change one number and high dosing stops working, which is exactly what the mutant prevention concentration literature has always said [5][6].

Mutation supply. At \(\mu = 10^{-8}\) the long arm ends in resistance 1000 times out of 1000, and at \(10^{-10}\) it cures 78.80%. A factor of a hundred in a parameter nobody can measure precisely for a given patient moves the answer across its whole range.

What the model does not contain

One compartment, so no tissue penetration and no site where the drug arrives weak. No immune system beyond the crude constant term of §10. One resistance mutation with one fixed cost, so none of the stepwise ratcheting through intermediate MIC levels that real fluoroquinolone resistance shows [1][7], and no compensatory mutations eroding the cost over time [9]. No back mutation and no horizontal transfer, by plasmid or otherwise. A well-mixed population, so no biofilm and no spatial refuge. No persisters and no tolerance. Identical pharmacokinetics in every replicate, where real between-patient variation in clearance is large. Adherence that is independent per dose, where real non-adherence comes in runs. And one patient at a time, with nothing at all about transmission of resistant strains between people, which is the part of this problem that actually matters for public health [10].

Any one of those could reverse a result here. The one we would look at first is the immune term, for the reason given in §10, and the second is the bursty structure of real non-adherence, because §9 says that the shape of the gaps is what does the damage and we modelled the gaps in the least realistic way available.

Reproducing this

Everything in this article comes from one file and one command. Python 3.12 and numpy are the only requirements.

python resistance-evolution.py > resistance-evolution-output.txt

Expect six to nine minutes. Ours took 498.0 seconds on a laptop under Python 3.12.3 and numpy 2.4.2, and produced 178,002 simulated treatment courses. The master seed is 20250915, hard-coded at the top of the file, and every cell draws an independent stream from it through numpy's SeedSequence.spawn, so the output is deterministic and you should get our numbers to the last digit. If you do not, we would like to know. The interactive model runs a smaller version of the same simulation in your browser and reproduces the headline numbers on its default settings.

A last word, since the subject invites misreading. This article is about a program we wrote. It compared schedules inside that program and found that peak height and the shape of the gaps decided the outcome while total drug and course length mostly did not. Whether any of that is true of a person with an infection is a question for clinical trials, some of which have been asking it [17][18], and for the doctor treating that person. Do not change how you take a medicine because of a simulation written by a school club.

References