# A model-based approach to characterize enzyme-mediated response to antibiotic treatments: going beyond the SIR classification

Virgile Andreani, Lingchong You, Philippe Glaser, Gregory Batt

► To cite this version:

Virgile Andreani, Lingchong You, Philippe Glaser, Gregory Batt. A model-based approach to characterize enzyme-mediated response to antibiotic treatments: going beyond the SIR classification. 2021. [hal-03544031](https://hal HAL Id: hal-03544031)

https://inria.hal.science/hal-03544031

Preprint submitted on 28 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A model-based approach to characterize enzyme-mediated response to antibiotic treatments: going beyond the SIR classification

2

3

Virgile Andreani¹,²*, Lingchong You³, Philippe Glaser⁴ and Gregory Batt¹,²*

4

1. Inria Paris, 2 rue Simone Iff, 75012 Paris, France

5

2. Institut Pasteur, 28 rue du Docteur Roux, 75015 Paris, France

6

3. Department of Biomedical Engineering and Center for Genomic and Computational Biology, Duke University, Durham, NC, USA

7

8

4. Unité EERA, CNRS UMR 3525, Institut Pasteur, APHP, Université Paris-Saclay, 28 rue du Docteur Roux, 75015 Paris, France

9

10

* Corresponding authors: virgile.andreani@inria.fr, gregory.batt@inria.fr

11

## Abstract

12

13

14

15

16

17

18

19

20

21

22

23

To design appropriate treatments, one must be able to characterize accurately the response of bacteria to antibiotics. When exposed to β-lactam treatments, bacteria can be resistant and/or tolerant, and populations can exhibit resilience. Disentangling these phenomena is challenging and no consolidated understanding has been proposed so far. Because these responses involve processes happening at several levels, including the molecular level (e.g. antibiotic degradation), the cell physiology level (filamentation) and the population level (release of β-lactamases into the environment), quantitative modelling approaches are needed. Here, we propose a model of bacterial response to β-lactam treatments that accounts for bacterial resistance, tolerance, and population resilience. Our model can be calibrated solely based on optical density readouts, can predict the inoculum effect, and leads to a mechanistically relevant classification of bacterial response to treatments that goes beyond the classical susceptible / intermediate / resistant classification. Filamentation-mediated tolerance and collective enzyme-mediated antibiotic degradation are essential model features to explain the complex observed response of cell populations to antibiotic treatments.

24

## Introduction

25

26

27

28

29

30

31

32

33

34

35

Antibiotic resistance is one of the major health concerns of our time. To enable appropriate treatment decisions for specific clinical isolates, one needs to characterize them as precisely as possible. The standard susceptible / intermediate / resistant classification (SIR) relies on the measure under standard conditions of a single indicator of susceptibility of a strain to an antibiotic, the minimal inhibitory concentration (MIC). However, the response of a cell population to an antibiotic treatment is multifactorial and can be described in terms of resistance (the ability of a cell to *grow and divide* in presence of antibiotics), tolerance (the ability of a cell to *survive* in presence of antibiotics), or resilience (the ability of a cell *population* to *recover* after a treatment). This raises the question whether one can extract more information about the susceptibility of cells to an antibiotic from the whole growth curve of a cell population (Figure 1A). This is all the more relevant given that optical density (OD) growth curves can be easily accessible for diagnostic purpose in clinical environments.

36

37

Although the fine details of the mechanisms are still being investigated¹⁻⁵, it is known that the molecular mechanism of action of β-lactam antibiotics is the inactivation of penicillin binding proteins (PBPs), which
disables the capacity of a cell to synthesize or repair its cell wall. Earlier studies have also shown that these
molecular actions can lead to dramatic cell shape changes in bacteria under antibiotic treatment. Many
β-lactams cause filamentation, but some can also cause the apparition of spheroplasts, ovoid or bulging
cells⁶⁻⁸. These processes delay lysis and hence are a form of tolerance⁹. Research has also been done on the
response of a cell population to β-lactam treatments, highlighting the role of enzyme-mediated resilience (also
called collective antibiotic tolerance (CAT)) either through OD measurements or cell counting¹⁰⁻¹². However,
even if escape mechanisms are now partly understood at the molecular, cellular, and population levels, no
consolidated vision of their combined effects has been proposed so far to quantitatively explain the response
dynamics of cell populations treated with β-lactams. We attribute this to the complexity of the phenomena
involved, but also to the intrinsic experimental or conceptual problems coming with the two main observables
of this system: cell countings requires constant and laborious interventions, and the interpretation of OD
measurements is notoriously difficult in presence of filamentation or swelling.

In this study, we tackle the challenging problem of the precise characterization of strains through OD
measurements only, to enable the approach to be implementable in clinical settings that are often limited to
this kind of measurements. Thus, we focused on the use of a plater reader to acquire OD response curves of
various *E. coli* clinical isolates to β-lactam treatments in micro-plate wells. We observed a wide variety of
response dynamics. Yet, limiting the analysis to this data poses an informational challenge notably because of
the decorrelation between live cell numbers and OD readouts in presence of filamentation, and of the difficulty
to access other quantities of interest in the system like the evolution of the concentrations of chemical species
in the culture medium. For this reason, we assembled the knowledge of the field on bacterial response to
β-lactams into a mathematical model encompassing the molecular, cell physiology and population levels,
general enough to capture the dynamics of all considered isolates, and simple enough to be identifiable only
from OD data. This dedicated modelling effort resulted in a model able to correctly predict well known
observations such as the inoculum effect, which describes how the MIC increases with the initial number of
cells. Fitting this model to the data we collected on clinical isolates allowed us to extract parameter values
and uncertainties specific to each individual isolate. This then led to the emergence of a natural classification
of the different isolates in terms of resistance and tolerance capacities. Because this classification makes
mechanistic sense, it goes beyond the SIR system. The full approach is summarized in Figure 1B.

Results

Capturing the dynamics of optical density, cell number and length distribution
in a mechanistic model

The challenge of the approach outlined in Figure 1B consists in characterizing many different processes
from incomplete but longitudinal data collected from cell cultures treated with antibiotics. To assess the
feasibility of this problem, we constructed a mathematical model of the response of a bacterial population to
β-lactam treatments targetting cell elongation, involving all relevant processes at the molecular, single-cell, and
population levels. Notably, the model incorporates also dynamical processes that are not readily measurable
in a clinical context and can therefore be used to reconstruct hidden dynamics from experimental readouts.

Modelling cell length and filamentation is needed in order to explain tolerance to β-lactam treatments⁹.
However, simply considering that all the cells have the same length, which evolves with time because of
filamentation, cannot explain the complex dynamics of cell lysis. Indeed, only the longest cells die, which
creates a non-trivial feedback on the average length of surviving cells. Hence, accounting for the heterogeneity
in cell lengths in the surviving cell population is required. Therefore, we developed a partial differential
equation (PDE) model of the temporal evolution of the distribution of cell lengths within the population,
and described with ordinary differential equations (ODEs) other processes such as the active degradation
of the antibiotic by β-lactamases released in the medium upon cell lysis. The evolution of the cell length
distribution is represented by a process of growth-fragmentation that relies on the assumption that cells
elongate exponentially and divide with a rate depending on their length, into possibly more than two fragments,
A

Problem

<figure><img src="image_2.png" /></figure>

B

Approach

<figure><img src="image_3.png" /><figcaption>**Figure 1: Characterizing antibiotic response using OD growth curves.** **A.** Exploiting all the information available in OD growth curves has the potential to be more informative. Top: Minimal inhibitory concentrations (MICs) are measured with a unique time point at 18 hours. Isolates A and B have a high MIC, whereas isolate C has a low MIC. Bottom: Following the optical density (OD) of the cultures during the experiment allows one to unveil that strains A and B behave differently despite having the same MIC. One can qualify A as highly resistant, and B as poorly resistant but tolerant and resilient. Strain C is neither resistant nor tolerant. **B.** Our approach. A 96-well plate is inoculated with the isolates to be tested, and a range of antibiotic concentrations is applied. The OD is read by an automated plate reader and processed automatically. The data is then fed into the model, for which numerical parameters are estimated, along with their uncertainties, for each isolate. On the basis of these parameters, a machine learning approach can be used to classify the strains into clusters.</figcaption></figure>
which is notably the case for filamented cells after removal of the antibiotic¹³. Because some β-lactams like cefotaxime can inhibit PBP3, involved in the septum formation¹⁴,¹⁵, the division rate is assumed to be a decreasing function of the antibiotic concentration in the medium. However, the elongation rate is only a function of the nutrient concentration, and remains unchanged by the antibiotic, as is the case for most β-lactams¹⁰.

It is known that for a given β-lactam concentration, the time to cell lysis is inversely proportional to the growth rate¹⁶. This suggests that the death mechanism can be easily modelled with a critical cell length, reachable when cell division is blocked. Cell filaments longer than this threshold experience a significant death rate (Figure 2A). This makes cell lysis a direct function not of the antibiotic itself but of cell length, that drastically increases because of the antibiotic. The antibiotic can however, at high doses, reduce this critical length, through the inhibition of other PBPs such as the PBP1s which have a fundamental role in the repair of wall defects⁵.

All the bacteria in this study express β-lactamases, enzymes released in the culture medium upon cell lysis, with the ability to cleave the antibiotic molecules¹⁷. Although strains can express several β-lactamases with different efficiencies, the model aggregates them into a single average one. The concentration of antibiotics, as well as of β-lactamases, are tracked by means of ODEs. The model also accounts for the excess OD contributed by fragments of lysed cells (the dead biomass). The equations of the corresponding PDE model are given in the Methods section and are described in more detail in the Supplementary Text A and B.

To make the model computationally tractable, we eliminated the explicit representation of the cell length distribution $n(\ell, t)$ to keep only its first moments: number of cells $N$ and average length $L$. The ODEs for these quantities involve partial moments of the cell length distribution, which we managed to express only in terms of $N$ and $L$ through careful approximations (see Supplementary Text A for details).

Automating a microplate reader with the custom driver `plateRider`, we generated experimental data using a protocol designed to exhibit a strict Monod growth in absence of treatment, to limit evaporation and to correct for non-linearities of OD measurements (see Methods and Supplementary Text C, D, E and F). The model was fitted to OD and live cell counting data as described in Methods. Figure 2B shows the experimental data for a reference strain subjected to different initial concentrations of cefotaxime, and the model simulations fitted to the data. We show on this figure that our model offers a very good agreement with the experimental data and additional insights about processes that are not easily measurable such as the average cell length or the chemical composition of the culture medium.

Furthermore and importantly, the agreement remains good whether we fit it on OD and cell number, or OD only and make it predict the cell number (see also Supplementary Text G for another isolate). This shows that there is enough information in the OD growth curves to lift the veil on important hidden processes such as the evolution of the live cell number. Lastly, it also qualitatively captures the evolution of the cell lengths, observed with a microscope (Supplementary Text H). To further challenge the predictive power of our model, we tested its capacity to predict quantitatively the inoculum effect.

## Prediction of the inoculum effect

Active antibiotic degradation by β-lactamases released in the culture medium upon cell lysis is known to cause a decrease of the observed antibiotic efficacy as the initial cell density increases, a phenomenon known as the inoculum effect¹⁸,¹⁹ (Figure 3A). This might lead to the regrowth of the population, then showing resilience, also known as collective antibiotic tolerance. The quantitative characterization of this effect has been the subject of several studies¹¹,²⁰, as its detailed understanding would be crucial in the design of optimal treatments²¹, or treatments adapted to the estimated cell density at the site of infection²².

A reference strain (ANSES 32139) was grown in presence of a range of antibiotic concentrations, starting from several cell densities. As shown in Figure 3B (left), the inoculum effect emerges indeed from simulations of the model, even though the model was in this case only fitted on two of the eight seeding densities and is free to extrapolate on the others. Obtaining these predictions is far from trivial as it requires not only
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2: A multiscale model of β-lactam treatments on cell populations. A.** Graphical representation of the model. *Cell response:* Cells elongate exponentially at a rate *g* that depends on nutrients but is not affected by the antibiotic. Above a concentration *k*₁, the antibiotic inhibits PBP3, which deprives the cells of their capacity to divide. Cells that cannot divide filament until they reach a critical length *L*<sub>max</sub> from where they experience a cell death with rate γ. Higher concentrations of antibiotics, around *k*₂, inhibit the PBP1s which fragilises the cells and reduces the critical length, such that cells die earlier. On the right are pictured schematic representations of the typical response of the cell culture to an initial antibiotic treatment. *Population response:* When the cell population is attacked by the β-lactam treatment (red dots), cells elongate and eventually lyse, releasing β-lactamase in the environment (blue circles) that degrade the antibiotic thus protecting the surviving cells. **B.** Comparison of experimental data and the model fit. The isolate ANSES 32139 is treated with different concentrations of cefotaxime (CTX), and optical density and live cell numbers (dots) are measured as described in the Methods section. For cell numbers we also show 95% confidence intervals on the CFU inverse problem estimation (see Methods). Data for OD has been subsampled for clarity. We show in solid lines the output of a model fitted on OD only, and in dashed lines a model simultaneously fitted on OD and cell number. The model also outputs predictions for the average cell length, and the concentration of antibiotic and other chemical species in the culture medium.</figcaption></figure>
132

133

134

capturing accurately the mechanism of clearance of the antibiotic by the β-lactamases and the evolution of the number of surviving cells to predict the regrowth of the population, but also extrapolating this information to conditions that the model never saw.

135

136

137

To further characterize the effect, the experimental, and predicted MICs were plotted together in Figure 3B (right). In accordance with recent experimental as well as theoretical results¹¹, we found that at low cell densities, the measured MIC is an exponential function of the initial cell density.

A

Principle of the inoculum effect

<figure><img src="image_6.png" /></figure>

B

Predicting treatment effects for different seeding densities

<figure><img src="image_7.png" /><figcaption>**Figure 3: Predicting quantitatively the inoculum effect.** A. The effect of antibiotic depends on the initial cell density. For low cell densities, the released β-lactamases are insufficient to clear the antibiotic of the medium before the lysis of all of the cells. For high cell densities, enough “sacrificial cells” can die to release a quantity of β-lactamase able to clear the antibiotic before all cells lyse. The remaining cells can then eventually divide and reconstitute the population, showing resilience. B. On the left, OD growth curves measured for different inocula are predicted by the model, after a fit on only two intermediate inocula. Both the crash and the regrowth, yet different for each inoculum / antibiotic concentration pair, are accurately predicted. On the right, the evolution of the MIC as function of the initial cell number is also predicted. Vertical bars represent the lowest (highest) concentrations of antibiotics tested that produced a cell OD lower (higher) than 0.1 after 20 hours.</figcaption></figure>

138

Capturing the diversity of responses of a selection of clinical isolates

139

140

141

142

Growth curve experiments were carried out on a collection of 9 *E. coli* clinical isolates selected to express a representative panel of different families of β-lactamases, including carbapenemases. The isolates, all belonging to the sequence type (ST) 410, have been chosen with sensitivities to cefotaxime ranging from fully susceptible to highly resistant. Moreover, several of these isolates contain mutations in PBP3 (*ftsI*) or porins
(ompC & ompF) contributing to a decreased antibiotic susceptibility²³ (Figure 4A and Supplementary Text I).

The same parameter fitting procedure was applied to each isolate shown in Figure 4B. We obtained very good fits to the data with our model, a result which we found challenging to reproduce with any other standard ODE model that we considered. We concluded that capturing the heterogeneity of cell lengths, either explicitly using a PDE model or implicitly with the ODE model derived from it, is essential. Moreover, our model quantitatively captures not only the OD curves but also the evolution of the number of cells, as shown for two of the nine isolates in Figure 2B and Supplementary Text G.

Interpreting resistance-conferring mutations as discrete events significantly changing the value of the associated parameters, we looked for signs of the mutations of *ftsI*, *ompC* and *ompF* in the parameter sets corresponding to the isolates where they appear. We found indeed that the highest values of the parameter controlling the sensitivity of PBP3 to the antibiotic were attributed to the isolates with mutations on *ftsI*, in the expected order from the effectiveness of the mutations (Supplementary Text J).

Furthermore, using literature data, we estimated the effectiveness of the CTX-M-1, CTX-M-15, CTX-M-55, TEM-1, NBM-5, CMY-2, CMY-22, and OXA-181 β-lactamases on cefotaxime, and found that isolates with several β-lactamases were indeed associated to models with high values for the parameters that control the β-lactamase efficiency (Supplementary Text K).

## Classification of isolates in meaningful categories by identifiability analysis

Besides their intrinsic value describing the strain, the nine parameter sets obtained by parameter estimation from the nine isolates can also be considered in relation to one another. Similar behaviours in different isolates should manifest in similar parameter values among these strains. To avoid the visualization of data in a 17-dimensional space, we used a t-distributed stochastic neighbour embedding (t-SNE), a statistical method to visualize high-dimensional data on low-dimensional maps. t-SNE respects proximity between points: two points close in the high-dimensional space should be mapped to points close in the low-dimensional space.

However, direct clustering of isolates based on parameter values fails to robustly identify clusters (Supplementary Text L). This is attributable to the presence of unidentifiabilities in the model. Even though we verified that the model is in principle fully identifiable, it can always happen that the range of antibiotics used experimentally is not sufficient to induce a strain to fully exhibit all the possible expected behaviours. Hence, for highly resistant or for susceptible isolates, the available data is not sufficient to allow the extraction of a well-constrained value for all model parameters. In these cases, these parameters are unidentifiable in practice, and their fitted values suffer from large variabilities. Therefore, two related isolates might not look similar in terms of parameter values because of the unidentifiable part of the model.

However, we realized that this identifiability challenge can be exploited. Indeed, similar isolates might be more similar in parameter *identifiabilities* than in parameter *values*. Parameter fitting algorithms can also return, additionally to the parameter set maximizing the likelihood of the data, an estimation of the curvature matrix of the likelihood at the optimal point. In first approximation, this curvature matrix can be used to construct confidence intervals around the optimal parameter set.

To first test this approach *in silico*, we generated a collection of 515 synthetic clinical isolates by random sampling of parameter values within biologically plausible limits. We used these parameters to generate data with simulated experimental noise and used our calibration approach on these data to obtain fitted parameter values and uncertainties. As seen for real experimental data, clustering the synthetic isolates on the basis of parameter values does not lead to conclusive results (Supplementary Text L). In contrast, clustering the synthetic isolates on the basis of these confidence intervals yields several groups that are phenotypically homogeneous, as shown in Figure 4C (left). This classification appears mostly orthogonal to SIR, several groups being constituted of a mixture of S, I and R isolates.

Treating the collection of 9 clinical isolates in the same way, we observed the emergence of three clusters.
A

Genetic features of selected clinical isolates

<table><thead><tr><th></th><th>Strain</th><th>CTX-M-1</th><th>CTX-M-15</th><th>CTX-M-55</th><th>TEM-1</th><th>NBM-5</th><th>CMY-2</th><th>CMY-42</th><th>OXA-1</th><th>OXA-181</th><th>ftsI</th><th>ompC</th><th>ompF</th></tr></thead><tbody><tr><td>1</td><td>ANSES 28668</td><td>x</td><td></td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>ANSES 32139</td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td><td>x</td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>ANSES 29401</td><td></td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>4</td><td>ANSES 30599</td><td></td><td></td><td></td><td>x</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>5</td><td>CNR 49A5</td><td></td><td>x</td><td></td><td>x</td><td>x</td><td></td><td>x</td><td>x</td><td></td><td></td><td></td><td>YRIK</td></tr><tr><td>6</td><td>CNR 94G8</td><td></td><td>x</td><td></td><td></td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td>YRIK</td></tr><tr><td>7</td><td>CNR 82A2</td><td></td><td></td><td>x</td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td>x</td><td>YRIN 349-532</td><td>R195L -46; C-&gt;T</td></tr><tr><td>8</td><td>CNR 84G4</td><td></td><td></td><td>x</td><td>x</td><td></td><td>x</td><td></td><td>x</td><td></td><td>x</td><td>YRIN 349-532</td><td>R195L -46; C-&gt;T</td></tr><tr><td>9</td><td>CNR 92B5</td><td></td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td>x</td><td></td><td>YRIN 349-532</td><td>R195L -46; C-&gt;T</td></tr></tbody></table>

B

Model fitting

<figure><img src="image_8.png" /></figure>

C

Strain classification based on fitted parameters

<figure><img src="image_9.png" /><figcaption>**Figure 4: Characterizing clinical isolates.** A. Each clinical isolate considered in this study codes for one or more β-lactamases, including carbapenemases. Five of them also show mutations on genes coding for PBP3 or porins, decreasing the antibiotic susceptibility of the cell. B. Experimental data has been collected for the 9 clinical isolates with 12 initial antibiotic concentrations each (left) and parameter fitting was applied to OD measurements (right). C. t-SNE clustering based on the confidence intervals of fitted parameters. Left: 515 random parameter sets were used to generate synthetic experiments, including simulated measurement noise. We fitted the model on these experiments and classified the synthetic isolates according to the practical identifiability of these fitted parameters. Several groups emerged, that are phenotypically homogeneous. Groups that do not exhibit growth limitations are shown more transparent as they correspond to implausible biological conditions. Right: the 9 clinical isolates were classified in the same way, creating three of the phenotypical groups observed in the synthetic approach. The two isolates labelled N are fitted on both OD data and number of cells, the rest is only fitted on OD data. On both sides, one representative plot is provided for each group.</figcaption></figure>
Importantly, each cluster can be associated to a specific way to respond to treatment. More specifically, we obtain one cluster for strains that do not show significant filamentation and are not affected by the antibiotic, one cluster for strains that are impacted by the antibiotic with filamentation and massive cell death but that are tolerant enough so that a small sub-population of survivors recovers from intermediate treatments, and one cluster for isolates that filament, die earlier, and cannot recover even from intermediate treatments. That is, using a non-supervised clustering approach, we observed the *emergence* of a classification of isolate responses into three mechanistically relevant classes: susceptible, tolerant, and resistant (STR). One should stress that this classification is different from the SIR system: isolates classified as resistant based on their MIC values are distinguished based on their actual response to treatments into tolerant and *bona fide* resistant isolates. Interestingly, the same classes were present among the groups observed in the synthetic study: resistant isolates of the collection were found in the H cluster, tolerant isolates were found in the E cluster, and susceptible isolates in the C cluster. Yet, not all behaviours observed in our collection of synthetic isolates have been observed in a real isolate. This raises the question whether these behaviours are biologically impossible or are simply not present in our collection.

## Discussion

With medicine increasingly turning towards personalization, the rapid identification of the properties of an isolate causing an infection represents a concrete and useful information for the clinician treating it. The current standard is a 3-group system: SIR, standing for susceptible, intermediate and resistant bacteria. The inclusion criteria into these groups consist of comparing the MIC of an isolate to thresholds, specific to cell types and antibiotics, updated every year by health authorities. However, describing the response of an isolate to treatments by a unique number leads necessarily to an extremely coarse representation. Moreover, the distinction between the different groups relies on rather arbitrary threshold values with no mechanistic interpretation.

In this study, we presented a model-based approach to classify *E. coli* clinical isolates exhibiting enzyme-mediated antibiotic resistance, tolerance and population resilience. The model, based on a growth-fragmentation equation and on a limited number of simple biological hypotheses, represents a comprehensive understanding of phenomena happening at the molecular, cell physiology, and population levels, as well as their interactions. Notably, it gives a central role to antibiotic-mediated filamentation as a way for cells to gain time until degradation of the antibiotic.

To our knowledge, our model is the first model of bacterial resistance, tolerance and resilience to β-lactams able to predict simultaneously the number of cells and optical density of a cell culture subject to filamentation and submitted to complex antibiotic treatments. It solves an observational and modelling challenge never tackled head-on by the community. Previous studies indeed avoided the filamentation issue by focusing exclusively on cell number measurements²⁴,²⁵, or when OD was considered, the pre-crash portion of the OD growth curve had to be ignored if the model did not account for filamentation¹²,²⁶. Ignoring the pre-crash portion of the growth curve prevents notably to quantify the effect of the antibiotics on PBP1s.

In addition, we demonstrated that our model is able to predict quantitatively the inoculum effect over a wide range of inocula, relying on a characterization done only on a small range. Our modelling of the inoculum effect emerges from a model of cell response which was not specifically developed for the inoculum effect, and is thus not an *ad hoc* phenomenological model. Our model also leads to parameter estimates consistent with the known genotype of the isolates.

Importantly, the proposed approach is powerful enough to only require OD data for model calibration, an observable easily accessible but notoriously hard to exploit in presence of filamentation²⁷. We used this feature to obtain information on parameter values and sensitivities for a selection of diverse clinical isolates, and used a non-supervised clustering approach to classify their behaviours. We observed the emergence of three distinct groups with mechanistic interpretations, leading to the proposition of a susceptible, tolerant and *bona fide* resistant strain classification, a classification based on a significantly richer information than the SIR
system. This result confirms the presence of bacterial identity and traits in their growth dynamics²⁸, and by comparison with a model-free visual or automated classification on raw OD curves, adds an explanatory layer that consists of a limited number of biologically-sensible model parameters, explaining escape mechanisms.

One aspect of the response of bacterial populations to antibiotics not included in the model is persistence, an extreme case of tolerance used to describe a subpopulation of cells with a distinctively different phenotype than the rest of the population. Persistent cells do not grow, which is why they are not affected by length-induced lysis. The current understanding of persistors is that they can survive during a long antibiotic exposure, and switch stochastically to the normal phenotype to reconstitute the population after the drug has been cleared. The fact that persistors are not needed in our model to explain the regrowth of the population shows that sufficiently tolerant isolates can also exhibit the ability to “survive exposure to high concentrations of an antibiotic”, a trait usually associated with persistence²⁹.

Our model was developed for entero-bacteria submitted to treatments with cefotaxime. This antibiotic primarily targets PBPs, proteins involved in cell division, elongation, and maintenance of cell shape. However, different β-lactams might have different affinities to different PBPs, and inhibit them in different orders, leading to different physiological responses of cells. For example, while ampicillin and cefotaxime, preferentially binding to PBP3, trigger a marked filamentation response, mecillinam, with a greater affinity to PBP2, is known to cause the apparition of ovoid cells. Even though our model assumes that filamentation is the prime mechanism of tolerance, the method is in principle applicable more widely. Thus, we anticipate that only slight adaptations of the model will be needed to adapt it to other β-lactam antibiotics.

Our model captures not only the cell response to treatment (pharmacodynamics) but also the evolution of the cells' environment, that is, of the culture medium in the well of a microplate (elementary pharmacokinetics). This work could therefore be extended in two directions. Firstly, it could be interesting to grow cells in low-volume parallelized bioreactors to be able to cultivate cells over longer durations and at higher densities, but also to observe the impact of persistors and of mutations, and relate observations with the model parameters and the STR classification of the isolates. Secondly, by completing our model with a proper pharmacokinetics component, we would obtain an excellent starting point to quantitatively investigate the effects of treatments *in vivo*, in which our model could notably help test whether the local release of β-lactams at the site of infection can have a real impact on treatment outcome. There is indeed a known correlation between the presence of an *in vitro* inoculum effect and failures of *in vivo* treatments²².

# Materials and Methods

## Strains and antibiotics

The *E. coli* sequence type 410 isolates used in this study are described in brief in Figure 4A and in more details in Supplementary Text I. They have already been the object of previous studies²³. Cefotaxime was purchased from Sigma-Aldrich and dilutions were made taking into account the purity specified by the vendor. Antibiotics and stock solutions were kept at -20°C and renewed regularly to limit degradation.

## Growth conditions

Unless specified otherwise, all pre-cultures, cultures and experiments were performed in 0.1% glucose M9 liquid medium (1 g/L glucose, 6.78 g/L Na₂HPO₄, 3 g/L KH₂PO₄, 1 g/L NH₄Cl, 0.5 g/L NaCl, 0.24 g/L MgSO₄, 0.01 g/L CaCl₂). The low glucose concentration is meant to create the conditions of a carbon-related growth arrest, a situation that lends itself better to mathematical modelling (Monod growth, see Supplementary Text D). The growth arrest happens in this medium between 0.2 and 0.3 OD₆₀₀.
## Growth curves acquisition

**Overnight.** A single bacterial colony was picked from an agar plate and incubated overnight at 37°C and 200 rpm. **Preculture.** The overnight was centrifuged, resuspended in fresh M9 and diluted to 0.05 OD$_{600}$ before a new incubation of 3 to 4 hours under the same conditions, aiming to catch the cells in exponential phase for the beginning of the experiment. **Experiment preparation.** Cells from the preculture were centrifuged, resuspended in fresh M9 and diluted to 0.05 OD$_{600}$. Each well of a 96-well plate with transparent flat bottoms was filled with 196 µL of M9, 2 µL of the cell suspension, and 2 µL of a stock antibiotic dilution. Depending on the plate configuration and the number of replicates, the dilutions were made as much as possible in larger quantities in order to minimize errors related to small-volume pipetting. **Data collection.** A 2017 Tecan Spark multimode plate reader was used for all OD acquisitions, in combination with a custom software driver, plateRider, giving us extensive control on measurements (Supplementary Text C). The plates used were Corning 3370 flat bottom microtiter plates. The 96-well plate was enclosed in a “humidity cassette” to limit evaporation. OD measurements were carried out in a loop consisting of measurement and incubation with shaking. The humidity cassette is only open during the measurements, following a schedule that was optimized to limit evaporation (Supplementary Text E). The incubation periods are carried out at 37°C and last 300 s, during which five incubation cycles happen (40 s of still incubation followed by 10 s of linear shaking (amplitude 1.8 mm and frequency 135 Hz) and 10 s of orbital shaking (amplitude 5 mm and frequency 20 Hz). At each data point, 20 OD measurements of the same well are taken in short succession, and the median OD is returned, to eliminate possible outlier measurements.

## Live cell number estimation

Live cell number was estimated by the CFU counting method, with a custom statistical analysis method allowing to estimate the expected number of cells and a 95% confidence interval on each value, with as few as one replicate. We used bayesian inference to determine that if we observe $O_1$ CFUs out of a sampled volume $V_0$ of a cell culture of volume $V_1$, then the probability that the initial number of cells is $N$ is given by

$$p(N_0 = N | N_1 = O_1) = \binom{N}{O_1} \left(\frac{V_1}{V_0}\right)^{O_1+1} \left(1 - \frac{V_1}{V_0}\right)^{N-O_1}$$

The formula can be extended to incorporate several replicas. This formula was used to determine the expected initial number of cells and 95% confidence intervals around this value (see Ref$^{30}$ for more details). This was instrumental in minimizing the number of plates spread and counted.

## Data analysis

**OD blanking.** The measured OD of a well is the sum of the OD of the cell culture, and of the OD of the well bottom, which slightly varies among wells of the same plate. For this reason, a constant OD was removed from each well data, such that the initial measured value matched the experimental initial cell dilution. **Removal of aberrant time points.** Occasionally, an OD reading fails and returns a value much larger than both the previous and the following ones of the same well. As this situation is highly biologically improbable, these outliers were removed before model fitting. **Correction of the non-linearity of the measurements.** The relation between measured and theoretical OD was studied and measurements were corrected for non-linearities (see Supplementary Text F).

## Models

To describe the temporal evolution of the distribution of cell lengths $n(l, t)$, we used a growth-fragmentation model. It relies on four main mechanisms: the cells continuously elongate with a rate $g$, they divide with a rate $f$ and when they divide, regardless of their length, they always split into cells of sizes $\frac{1}{2}$ to 1 (arbitrary
317 unit), meaning that filamented cells might split into more than 2 cells. Finally, they lyse when they reach the critical length $L_m$. The following PDE describes these phenomena:

$$ \frac{\partial n}{\partial t} + g \ell \frac{\partial n}{\partial \ell} + g n = \begin{cases} 0 & \text{for } 0 \le \ell < \frac{1}{2} \\ f \sum_{i=1}^{\infty} 2^{2i} n(2^i \ell, t) & \text{for } \frac{1}{2} \le \ell < 1 \\ -f n & \text{for } 1 \le \ell < L_m \\ -(f + \gamma) n & \text{for } L_m \le \ell \end{cases} $$

319 The PDE is initialized with the steady state cell length distribution of an exponentially growing cell population: $n(\ell, 0) = N(0) y_{\infty, \gamma=0}(\ell)$ with

$$ y_{\infty, \gamma=0}(\ell) = \begin{cases} \frac{f+g \ell^{f/g} - 2^{-f/g}}{f \ell^{2+f/g}} & \text{for } \frac{1}{2} \le \ell < 1 \\ \frac{f+g (1-2^{-f/g})}{f \ell^{2+f/g}} & \text{for } 1 \le \ell \end{cases} $$

321 The interaction with the antibiotics is taken into account through the dependency of these rates and threshold with the antibiotics concentration in the culture medium. While the elongation rate is independent on the antibiotic concentration and just follows Monod's law, the division rate $f$ and the critical length $L_m$ are decreasing functions of the antibiotic:

$$ g = \mu \frac{s}{K_s + s} \qquad f = \frac{\beta}{1 + \left(\frac{a}{k_1}\right)^{h_1}} \qquad L_m = L_{\min} + \frac{L_{\max} - L_{\min}}{1 + \frac{a}{k_2}} $$

325 Various concentrations are also tracked in the culture medium, such as the nutrients $s$, the antibiotics $a$, the β-lactamases $b$, and the dead biomass $c$ and $c_r$. Their dynamics are described with the following ODEs:

$$ \begin{aligned} \frac{ds}{dt} &= -\frac{g}{\lambda} \int_{0}^{\infty} \ell n \, d\ell \\ \frac{da}{dt} &= -k_{b} b a - d_{a} a \\ \frac{db}{dt} &= \gamma B_{\text{in}} \int_{L_{\text{m}}}^{\infty} \ell n \, d\ell - d_{b} b \\ \frac{dc}{dt} &= \gamma (1 - p_{c}) \int_{L_{\text{m}}}^{\infty} \ell n \, d\ell - d_{c} c \\ \frac{dc_{r}}{dt} &= \gamma p_{c} \int_{L_{\text{m}}}^{\infty} \ell n \, d\ell \end{aligned} $$

327 Finally, the optical density is proportional to the sum of the live and dead biomasses:

$$ OD(t) = \eta \left( \int_{0}^{\infty} \ell n \, d\ell + c(t) + c_{r}(t) \right) $$

328 This describes the PDE model. For computational efficiency, we derived and used for calibration an ODE model approximating the PDE model above (see Supplementary Text A for detailed explanations).
## Model fitting

330

Each parameter was restricted to a range of biologically plausible values. Depending on the role of the parameter and the size of the range, a change of variable was applied or not, to perform the search of this parameter in the linear space or in a logarithmic space. All 17 search ranges were then brought back to the interval [0, 10] and it is in this [0, 10]$^{17}$ space that the parameter search was performed.

331

332

333

334

The cost function computes the log-likelihood of the data assuming independent Gaussian measurement noise on each point. Special care was taken when searching on mixed datasets involving both OD measurements and cell counts, to balance the contribution of both data types. The integration of the ODE system was done with the `solve_ivp` method of `scipy`$^{31}$ with the method 'LSODA', and absolute and relative tolerances set to $10^{-9}$. For the synthetic experiments, an additive Gaussian noise was applied to the OD, with a standard deviation that depends on the OD as calibrated on the real data: $\sigma = 0.02 \cdot OD + 10^{-4}$.

335

336

337

338

339

340

Model fitting was performed with CMAES$^{32}$, with an initial $\sigma = 2$, except when a local search was desired. This was notably the case for the estimation of parameter uncertainties. Then, the `least_squares` solver of `scipy`$^{31}$ was used with the method 'trf'. This solver also returns curvature information around the optimum, which we interpreted to compute asymptotic confidence intervals for the fitted parameters.

341

342

343

344

## Data and code availability

345

All data, `plateRider` scripts and log files are available on zenodo: 10.5281/zenodo.5111026. All scripts needed for model calibration and figure plotting are available as a git repository at the address https://gitlab.inria.fr/InBio/Public/esbl-escape. The custom plate reader driver `plateRider` is an open source software available at https://gitlab.inria.fr/InBio/Public/platerider.

346

347

348

349

## Acknowledgements

350

The authors thank Arthur Carcano for help with the fitting of a previous version of the model, Andela Davidović for discussions about the model and help with image analysis, Sebastian Sosa Carrillo and Chetan Aditya for assistance in the laboratory, Jean-Yves Madek for giving access to the isolates ANSES and Thierry Naas to the isolates CNR, Jakob Ruess for comments on the manuscript, Adriana Chiarelli, Nicolas Cabanel and Isabelle Rosinski-Chupin for discussions.

351

352

353

354

355

This work was supported by ANR grants Seq2DiAg (ANR-20-PAMR-0010) and Anoruti (ANR-20-PAMR-0001), and by the H2020 Fet-Open COSY-BIO grant (grant agreement no. 516 766840).

356

357

## Author contributions

358

VA and GB designed the research. VA designed and performed the experiments, derived the mathematical models, developed the plate reader driver and analysed the data. LY discussed the project and provided expertise on escape mechanisms. PG advised the project, provided expertise on β-lactam resistance and for the choice of the isolates. GB supervised the research. VA wrote the article with significant input from all other authors.

359

360

361

362

363

## Conflict of interest

364

The authors declare no competing interests.

365
References

366

367

Chung, H. S. et al. Rapid β-lactam-induced lysis requires successful assembly of the cell division machinery. *Proceedings of the National Academy of Sciences* **106**, 21872–21877 (2009).

368

369

Yao, Z., Kahne, D. & Kishony, R. Distinct single-cell morphological dynamics under β-lactam antibiotics. *Molecular Cell* **48**, 705–712 (2012).

370

371

Cho, H., Uehara, T. & Bernhardt, T. G. β-lactam antibiotics induce a lethal malfunctioning of the bacterial cell wall synthesis machinery. *Cell* **159**, 1300–1311 (2014).

372

373

Wu, S. et al. β-lactam antibiotics stimulate biofilm formation in non-typeable haemophilus influenzae by up-regulating carbohydrate metabolism. *PLOS ONE* **9**, e99204 (2014).

374

375

Vigouroux, A. et al. Class-A penicillin binding proteins do not contribute to cell shape but repair cell-wall defects. *eLife* **9**, e51998 (2020).

376

377

Greenwood, D. & Eley, A. A turbidimetric study of the responses of selected strains of pseudomonas aeruginosa to eight antipseudomonal β-lactam antibiotics. *Journal of Infectious Diseases* **145**, 110–117 (1982).

378

379

Buijs, J., Dofferhoff, A. S. M., Mouton, J. W., Wagenvoort, J. H. T. & van der Meer, J. W. M. Concentration-dependency of β-lactam-induced filament formation in gram-negative bacteria. *Clinical Microbiology and Infection* **14**, 344–349 (2008).

380

381

Fredborg, M. et al. Automated image analysis for quantification of filamentous bacteria. *BMC Microbiology* **15**, 255 (2015).

382

383

Zahir, T. et al. Image-based dynamic phenotyping reveals genetic determinants of filamentation-mediated β-lactam tolerance. *Frontiers in Microbiology* **11**, 374 (2020).

384

385

Rolinson, G. N. Effect of β-lactam antibiotics on bacterial cell growth rate. *Journal of General Microbiology* **120**, 317–323 (1980).

386

387

Artemova, T., Gerardin, Y., Dudley, C., Vega, N. M. & Gore, J. Isolated cell behavior drives the evolution of antibiotic resistance. *Molecular Systems Biology* **11**, 822 (2015).

388

389

Meredith, H. R. et al. Applying ecological resistance and resilience to dissect bacterial antibiotic responses. *Science Advances* **4**, eaau1873 (2018).

390

391

Wehrens, M. et al. Size laws and division ring dynamics in filamentous escherichia coli cells. *Current Biology* **28**, 972–979.e5 (2018).

392

393

Spratt, B. G. Distinct penicillin binding proteins involved in the division, elongation, and shape of escherichia coli K12. *Proceedings of the National Academy of Sciences* **72**, 2999–3003 (1975).

394

395

Botta, G. A. & Park, J. T. Evidence for involvement of penicillin-binding protein 3 in murein synthesis during septation but not during cell elongation. *Journal of Bacteriology* **145**, 333–340 (1981).

396

397

Boman, H. G. & Eriksson, K. G. Penicillin induced lysis in escherichia coli. *Journal of General Microbiology* **31**, 339–352 (1963).

398

399

Livermore, D. M. β-lactamases: Quantity and resistance. *Clinical Microbiology and Infection* **3**, 4S10–4S19 (1997).

400

401

Vega, N. M. & Gore, J. Collective antibiotic resistance: Mechanisms and implications. *Current Opinion in Microbiology* **21**, 28–34 (2014).

402

403

Meredith, H. R., Srimani, J. K., Lee, A. J., Lopatkin, A. J. & You, L. Collective antibiotic tolerance: Mechanisms, dynamics and intervention. *Nature Chemical Biology* **11**, 182–188 (2015).

404

405

Salas, J. R., Jaberi-Douraki, M., Wen, X. & Volkova, V. V. Mathematical modeling of the 'inoculum effect': Six applicable models and the MIC advancement point concept. *FEMS microbiology letters* **367**, fnaa012 (2020).

406
407

21. Meredith, H. R., Lopatkin, A. J., Anderson, D. J. & You, L. Bacterial temporal dynamics enable optimal design of antibiotic treatment. *PLOS Computational Biology* 11, e1004201 (2015).

408

409

22. Lenhard, J. R. & Bulman, Z. P. Inoculum effect of β-lactam antibiotics. *Journal of Antimicrobial Chemotherapy* 74, 2825-2843 (2019).

410

411

23. Patiño-Navarrete, R. et al. Stepwise evolution and convergent recombination underlie the global dissemination of carbapenemase-producing escherichia coli. *Genome Medicine* 12, 10 (2020).

412

413

24. Nikolaou, M. & Tam, V. H. A new modeling approach to the effect of antimicrobial agents on heterogeneous microbial populations. *Journal of Mathematical Biology* 52, 154-182 (2006).

414

415

25. Clarelli, F. et al. Drug-target binding quantitatively predicts optimal antibiotic dose levels in quinolones. *PLOS Computational Biology* 16, e1008106 (2020).

416

417

26. Lee, A. J. et al. Robust, linear correlations between growth rates and β-lactam-mediated lysis rates. *Proceedings of the National Academy of Sciences* 115, 4069-4074 (2018).

418

419

27. Stevenson, K., McVey, A. F., Clark, I. B. N., Swain, P. S. & Pilizota, T. General calibration of microbial growth in microplate readers. *Scientific Reports* 6, 38828 (2016).

420

421

28. Zhang, C. et al. Temporal encoding of bacterial identity and traits in growth dynamics. *Proceedings of the National Academy of Sciences* 117, 20202-20210 (2020).

422

423

29. Brauner, A., Fridman, O., Gefen, O. & Balaban, N. Q. Distinguishing between resistance, tolerance and persistence to antibiotic treatment. *Nature Reviews Microbiology* 14, 320-330 (2016).

424

425

30. Andreani, V. Modelling and efficient characterization of enzyme-mediated response to antibiotic treatments. (Ecole polytechnique, 2020).

426

427

31. Virtanen, P. et al. SciPy 1.0: Fundamental algorithms for scientific computing in python. *Nature Methods* 17, 261-272 (2020).

428

429

32. Hansen, N. *The CMA evolution strategy: A comparing review.* (2006).

430