Whole-cell modeling of E. coli colonies enables quantification of single-cell heterogeneity in antibiotic responses

Antibiotic resistance poses mounting risks to human health, as current antibiotics are losing efficacy against increasingly resistant pathogenic bacteria. Of particular concern is the emergence of multidrug-resistant strains, which has been rapid among Gram-negative bacteria such as Escherichia coli. A large body of work has established that antibiotic resistance mechanisms depend on phenotypic heterogeneity, which may be mediated by stochastic expression of antibiotic resistance genes. The link between such molecular-level expression and the population levels that result is complex and multi-scale. Therefore, to better understand antibiotic resistance, what is needed are new mechanistic models that reflect single-cell phenotypic dynamics together with population-level heterogeneity, as an integrated whole. In this work, we sought to bridge single-cell and population-scale modeling by building upon our previous experience in “whole-cell” modeling, an approach which integrates mathematical and mechanistic descriptions of biological processes to recapitulate the experimentally observed behaviors of entire cells. To extend whole-cell modeling to the “whole-colony” scale, we embedded multiple instances of a whole-cell E. coli model within a model of a dynamic spatial environment, allowing us to run large, parallelized simulations on the cloud that contained all the molecular detail of the previous whole-cell model and many interactive effects of a colony growing in a shared environment. The resulting simulations were used to explore the response of E. coli to two antibiotics with different mechanisms of action, tetracycline and ampicillin, enabling us to identify sub-generationally-expressed genes, such as the beta-lactamase ampC, which contributed greatly to dramatic cellular differences in steady-state periplasmic ampicillin and was a significant factor in determining cell survival.


Introduction
This supplementary material includes a more detailed description of the E. coli model used for the simulations discussed in the main text. This model, which we call vivarium-ecoli, was a fork of wcEcoli -the Covert lab's E. coli Whole-Cell Modeling Project [23; 36]. It recreated the wcEcoli model using Vivarium -a framework for multiscale, compositional modeling [1]. This supplement describes the operation of the original model and how it was migrated to Vivarium (section 2), the new processes that turned the model into an agent which could be plugged into a spatial environment with other whole-cell models (section 3), and finally the implementation of antibiotic response mechanisms (section 4).

Single-cell E. coli model
The single-cell model used in this paper was based on a snapshot of our wcEcoli repository (https://github.com/CovertLab/wcEcoli) taken on May 20, 2021. Thus, we will begin by providing an overview of the wcEcoli model (2.1) before discussing the advantages offered by the new Vivarium-based model, including greater flexibility in process parameterization (2.2), data organization (2.3), process organization (2.4), simulation configuration (2.6), and simulation organization (2.7).

wcEcoli
The wcEcoli model integrated data from a wide array of databases and published reports to calculate a set of nearly 20,000 parameters for over 10,000 equations. These equations were divided into a set of modules called processes, each of which encapsulates a distinct cellular process (e.g. mRNA transcript initiation, mRNA transcript elongation, etc.). The model was initialized with counts for all molecules in the cell and used the equations contained within its processes to calculate how these counts change over time.
There were 12 main processes in vivarium-ecoli adapted from the wcEcoli model. A full description of most of these original processes and their operation can be found in the supplementary text for the original release [23]. Since these have already been described, we only list them in the current supplement. The only exception is the chromosome structure process, which was added after the original release and will be described in an upcoming publication. The code for all processes can be found under the ecoli.processes directory.
These are the processes migrated from the wcEcoli repository:

Simulation parameters
Most simulation parameters for wcEcoli were held in an object called sim data, which was generated from literature and database data by a pipeline called the parameter calculator or ParCa. Further details about the ParCa can be found in the supplementary material for Macklin et al. [23].
In wcEcoli, processes derived all their parameters from this centralized sim data object. By contrast, vivarium-ecoli parameterized processes using user-supplied dictionaries, providing significantly more flexibility in terms of how parameters are sourced. Indeed, in our antibiotic simulations, while most parameters were still drawn from sim data, others were loaded from JSON files (see 2.6) or calculated at run time (refer to ecoli.library.parameters module). To extract the relevant parameters for each process from sim data, we created a new LoadSimData class that can be found in the ecoli.library.sim data module.

Simulation state
In wcEcoli, the simulation state was stored in arrays that were designed to hold preset data types. By contrast, vivarium-ecoli uses dictionary-like internal states that impose no restrictions on the types of data that can be stored, accessed, and modified by processes.
Briefly, each node in the Vivarium simulation state was a Store object that held data. In a typical simulation step, processes received the current state of all connected stores, then returned a set of updates to be applied to those stores. The default file used to pre-populate stores at the start of a simulation (data/wcecoli t0.json) was created by saving the initial state of a wcEcoli simulation as a JSON file. Notably, vivarium-ecoli was not restricted to this single initial state file and could be initialized with any initial state in JSON format, including states saved mid-way through a cell cycle.

Save state
The ability to start a simulation from any initial state was new to vivarium-ecoli. To take advantage of this feature, our colony simulations can be configured to create a JSON containing the state of all cells and their shared environment at specific time points. This saved colony state could subsequently be loaded as the initial state of a future simulation and pick up from the saved point in simulated time. Note that certain internal states (e.g. the saved solution in the flux balance analysis solver used by metabolism, the state of the stochastic simulator used by protein complexation, etc.) could not be easily accessed, saved, or reloaded, preventing perfect reconstruction from saved states. In the simulations run for the present work, this imperfect reconstruction introduced slight differences between simulations that ran for the full 26002 seconds uninterrupted (e.g. the baseline simulations) and those that were initialized with a 11550-second saved state of a baseline simulation (e.g. the antibiotic simulations).

Pseudo-random number generator seeds
To ensure that simulations are reproducible, vivarium-ecoli allowed users to manually set the initial state and seed of simulations via JSON configuration files or command-line arguments (see 2.6 for more information on methods of simulation configuration). All simulations initialized with the same initial state and seed will generate the same results.
If the initial state contains only a single cell, the simulation seed was used to generate separate seeds for every pseudo-random number generator (PRNG) employed by that cell (most processes have their own PRNG). If the initial state contains multiple cells, the simulation seed was first converted to a series of seeds, one per cell, each of which was then used as described in the single-cell case. Refer to ecoli.library.sim data for the exact details of PRNG seeding.
Our model relied on PRNGs both for creating internal unique identifiers (IDs) and to model biological phenomena such as stochastic processes. Thus, to prevent the clashing of unique IDs after division, daughter cells were each assigned new simulation seeds by a PRNG in the mother cell. Similarly, when loading a saved state, an algorithm was used to avoid initializing the new simulation with the same seed as the simulation which generated the saved state (refer to the run simulation method of ecoli.composites.ecoli engine process).

Partitioning
Since the processes in our model share molecular resources, wcEcoli used a partitioning system to prevent "overdrafts" whereby processes collectively consume more resources than are present. Partitioning occurs at every time step as follows: 1. Each process reads the currently available molecule counts and requests counts based on what is available (calculate request). 2. The Allocator attempts to fulfill requests based on process priorities (see below). 3. Each process calculates a change in molecule counts based on the counts it was allocated (evolve state).
Most processes had equivalent neutral priority with the following exceptions: protein and RNA degradation had equally higher priorities, the two-component system process had lower priority, and transcription factor binding had the lowest priority. Chromosome structure and metabolism were the only processes that were not partitioned in the base model. Instead, they ran in succession after all other processes had finished updating the simulation state. Requests of higher-priority processes were always handled before the requests of lower-priority processes. When n processes of the same priority had requests r 1 , r 2 , ..., r n whose sum was greater than the unallocated count c u of a molecule, the Allocator allocated c i molecules to process i as follows, with random allocation of remainders caused by the floor function: In this way, even if all processes were to deplete their entire share of allocated molecules, the total count would remain non-negative, preventing overdrafts.
PartitionedProcess was the base class for Vivarium processes that were subject to molecular partitioning in wcEcoli. This class had calculate request and evolve state methods. Each partitioned process is wrapped with one Requester process and one Evolver process.
Requester processes called the calculate request of the wrapped process in the first step of the scheme above. Evolver processes called evolve state of their wrapped process in the third step.

Migration tests
To facilitate accurate migration of wcEcoli to the Vivarium framework, we developed a system for comparing each vivarium-ecoli process with its wcEcoli equivalent, ensuring that they function the same under a variety of conditions. This system included both unit tests and whole-model comparisons. Scripts containing migration tests for each process were placed in the migration directory. These tests were and continue to be run on every change made on the vivarium-ecoli repository.

EcoliSim
Vivarium-ecoli introduced a JSON-based interface for simulation configuration called EcoliSim, enabling the creation of customized configurations for different simulation runs. Using EcoliSim, simulations could be configured to add or remove processes from the base model, pass alternative parameters into given processes, and set the run-order of designated orderdependent processes (e.g. Metabolism), among other features. Configuration files could be written to "inherit" settings from other configuration files, or merged programmatically with other configurations using the EcoliSim interface. This allowed simulations to be built in a modular fashion, for example by combining several configurations that each add a few related processes and their associated parameters onto the base model.

EngineProcess
EngineProcess was a new Vivarium process created to improve the communication overhead of colony-scale simulations. This process served as a wrapper around an entire single-cell vivarium-ecoli model, allowing the Vivarium engine to assign an operating system (OS) process to each cell instead of each process within each cell. (Note that OS processes are distinct though conceptually similar to Vivarium processes.) The E. coli model had over a dozen processes that were required to communicate with one another for proper function (see 2.4). This communication was much faster when the processes could share memory instead of passing messages between OS processes. Thus, by packaging all the processes of individual cells within EngineProcess instances, we benefitted from fast communication between processes in the same cell without sacrificing the advantages of multiprocessing for multi-cell simulations. As of writing, there is already work underway to integrate this feature seamlessly into a future version of the Vivarium core software.

Cell division
Prior work proposed that E. coli may control cell size by elongating an approximately constant amount per cell cycle [5]. In accordance with this, each cell in our model initiated division upon reaching a dry mass equal to its initial dry mass plus a media-specific expected dry mass increase fitted by the parameter calculator (see 2.2). Noise was introduced by scaling this fitted dry mass increase by a factor randomly sampled from N (1, 0.1). In the rare case that the dry mass threshold was reached before the cell has accumulated at least two complete copies of its chromosome, division was delayed until chromosome replication had completed.
During division, all cell processes and stores were duplicated, resulting in one of each per daughter cell. Store values were divided as follows: • Bulk molecules: Each molecule went to each daughter cell with equal probability.
• Chromosomes: Each chromosome was assigned to a daughter cell with equal probability while ensuring that each daughter got at least one full chromosome. For simplicity, division happens instantaneously and the two resulting daughters are placed end-to-end in the spatial environment model.

Spatial environment model
The model represented the environment as a two-dimensional rectangular space. Agent locations were continuous coordinates (x, y) within this space, and each agent had an angle from the x-axis, length, width, thrust, torque, and mass. The environment was discretized into a lattice of sites, each with the same volume. The sites were larger than the agents to reflect how even strong concentration gradients yield approximately equal concentrations over the small distances bacteria span [24]. Each site in the lattice had a concentration for each molecule in the environment.
Importantly, the agents did not interact directly; instead, all of their interactions were mediated by changes in the environment, where agents were buffeted by physical forces, molecules diffused toward homogeneity, and the media could shift between environmental conditions.

Inputs and outputs
Inputs: Outputs: • Length • Outer membrane surface area • Inner membrane surface area • Periplasm volume • Cytoplasm volume

Model description
The physical shape of the cell was modeled by the newly introduced Shape process. During our simulations, the MassListener process periodically updated total cell volume by dividing current cell mass by the assumed constant density of 1.1 g/mL [2]. Shape then used the calculated volume to compute cell dimensions by assuming a capsule shape formed by a cylinder capped by hemispheres at each end. The process also assumed a constant width w, so cells grew exclusively by elongation [25]. Using this information, Shape computed the length l and outer surface area a o of the cell as follows: The process also used a parameter f p , the fraction of the cell's volume consumed by the periplasm, to calculate the volume of the periplasm v p , the volume of the cytoplasm v c , and the surface area of the inner membrane a i :

Parameters
Parameter Value Description f p 2 × 10 −1 Fraction of cell volume consumed by the periplasm [34]. w 1 µm Width of the cell (constant) [38]. Table SM2: Parameters for the shape process.

Multibody physics 3.2.1 Inputs and outputs
Inputs: Outputs:

Model description
The Multibody process was a wrapper around the physics engine pymunk (http://www. pymunk.org/), which can model individual cell agents as capsule-shaped rigid bodies that can move, grow, and collide. This engine was configured with elasticity (0.9) to simulate damped bacterial collisions, random jitter to model Brownian motion, and friction (0.9) to model cell-cell adhesion. For more information on the meaning of the elasticity and friction parameters, see the pymunk documentation. Multibody ran pymunk with a time step onetenth of its own two-second time step to simulate the movement of agents. It then updated each agent's location, angle, thrust, and torque. Upon division, the resulting daughter cells were placed end-to-end in the same orientation as the mother.
To simulate the low Reynolds environment bacteria experience [29], the process multiplied linear forces by d l and angular forces by d a every pymunk time step. Since d l , d a < 1, this multiplication mimicked the high drag of a low Reynolds number environment. Refer to Agmon et al. [1] for more details about the parameters chosen for this submodel.

Parameters
Parameter Value Description f j 1 × 10 −4 pN A random force applied to each agent to simulate Brownian motion. Manually tuned to recapitulate the behavior described in [30]. agent shape segment Assumed shape of each cell. Segments are cylinders capped with hemispheres.

Reaction diffusion 3.3.1 Inputs and outputs
Inputs: • Environmental concentrations • For each agent: -Location -Molecules to exchange with environment Outputs: • Environmental concentrations • For each agent: -External environment view -Molecules to exchange with environment

Model description
The ReactionDiffusion process simulated bounded two-dimensional fields of molecular concentrations. Each lattice site (x, y) held the local concentrations of any number of molecules, and diffusion simulated how they homogenized across local sites. At the beginning of each time step, before diffusion between sites was calculated, a set of user-defined, enzyme-catalyzed reactions could be simulated using Michaelis-Menten kinetics. The process assumed a Hill coefficient of n = 1 for all reactions. This was an appropriate assumption for the AmpC-catalyzed hydrolysis of ampicillin [26], which is described in 4.2.2.
Once reaction updates were applied, the process simulated diffusion by convolving (with reflection at the edges) the 2D Laplacian kernel over the lattice (C), multiplying by the diffusion constant D and the time step ∆t, and adding the result to the lattice: Each agent could uptake and secrete molecules at its position in the lattice. The model used the LocalField process to convert molecular exchanges between agents and the environment into concentration deltas that were applied at each agent's lattice site.

Antibiotic response model
This section describes all the new cellular processes that were added to model the effects of tetracycline and ampicillin exposure.

Membrane permeability 4.1.1 Inputs and outputs
Inputs: • Porin counts • Outer membrane surface area Outputs: • Permeabilities to tetracycline and ampicillin

Model description
The ease with which a molecule is able to diffuse across a membrane can be quantified as its permeability coefficient or permeability. The Permeability process computed permeabilities for tetracycline and ampicillin based on the abundance of OmpF, the porin primarily responsible for tetracycline [37] and ampicillin [21] trans-membrane diffusion.
Tetracycline can cross the outer membrane either through a porin, with permeability P outer,ompf,tet per unit of porin concentration in the membrane, or through the inner phospholipid bilayer with permeability P outer,bilayer,tet . With a porin count n ompf and outer membrane surface area of a o , the permeability for tetracycline crossing the outer membrane was calculated as: P outer,tet = n ompf a o P outer,ompf,tet + P outer,bilayer,tet Tetracycline can also cross the inner membrane, but it does not do so through porins. Since the permeability of the inner membrane to tetracycline is constant, it was not computed by the permeability process.
Ampicillin primarily enters cells through porins, not by traversing the phospholipid bilayer [21]. Its permeability through the outer membrane was computed as: P outer,amp = n ompf a o P outer,ompf,amp

Parameters
Parameter Value Description P outer,ompf,amp 6.63×10 −1 cm·µm 2 /s/porin Permeability of outer membrane to ampicillin [21] divided by the average simulated number of OmpF porins per square micron of the outer membrane. P outer,ompf,tet 2.35×10 −9 cm·µm 2 /s/porin Porin-attributable permeability of outer membrane to tetracycline [37] divided by the average simulated number of OmpF porins per square micron of the outer membrane. P outer,bilayer,tet 7 × 10 −8 cm/sec Permeability of outer membrane to tetracycline without porins [37]. Outputs: • Periplasmic and cytoplasmic antibiotic concentrations • Exchange of antibiotics with environment

Model description
The AntibioticTransportOdeint process simulates trans-membrane diffusion, export, and hydrolysis.

Tetracycline diffusion across the outer membrane
Under physiological conditions, tetracycline can form univalent cations by chelating magnesium ions that are readily available at the surface of the outer membrane (OM) [27; 37]. These magnesium-tetracycline chelates (Tc-Mg + ) preferentially cross the outer membrane by diffusing through the OmpF porin [37]. The rate of diffusion is thought to be dependent on the Donnan potential, which is generated by negatively charged, membrane-derived oligosaccharides in the periplasm that cannot diffuse through the OM [33].
To model the diffusion of this charged species across the outer membrane, we used the Goldman-Hodgkin-Katz (GHK) flux equation: where J S is the current density contributed by ion S (A/m 2 , positive means outward flow of positive charge), P S is the permeability of the outer membrane for ion S (m/s), [S] o is the external concentration of ion S (mol/L), [S] i is the internal concentration of ion S (mol/L), z s is the charge of ion S, E is the membrane potential (J/C), F is Faraday's constant (C/mol), R is the gas constant (J/K/mol), and T is the temperature (K).
The GHK flux equation assumes a constant electric field across the membrane and has previously been shown to fit biological diffusion data measured in vitro [10; 14; 16]. Prior mathematical analysis revealed that this "constant field" assumption holds only if the following conditions are met [22]: • The magnitude of the net charge density in the membrane must be small • The membrane must be thin Specifically, for membranes on the order of 100Å thick, including the outer membrane of E. coli [3; 11], the magnitude of the net charge density in the membrane must be lower than that produced by 10 −3 M of a univalent cation in solution. Since we added tetracycline at concentrations on the order of 10 −6 M, the constant field assumption was reasonable and the GHK flux equation could be safely applied.
To approximate the rate at which Tc-Mg + diffuses across the outer membrane, we performed a unit conversion on the GHK flux equation (Eq. 1) and rearranged to yield influx in the form of change in concentration of a molecule per unit time d[S] i dt : where v p is the periplasmic volume and a o is the outer membrane surface area.

Tetracycline diffusion across the inner membrane
Fick's first law of diffusion states that the rate of change in internal concentration d[S] i dt of a molecule with external concentration [S] o and permeability P S across a membrane with area A is: Since the inner membrane does not contain porins, tetracycline is believed to diffuse from the periplasm into the cytoplasm in an uncharged form [37]. This would mean that there are no electrical influences to consider, and Eq. 3 can be and was used instead of Eq. 2.

Ampicillin diffusion across the outer membrane
Ampicillin is an uncharged molecule, so its diffusion across the outer membrane was modeled using Fick's law (Eq. 3). Ampicillin diffusion across the inner membrane is minimal [21] and was not considered in our model.

Antibiotic export
Tetracycline and beta-lactams are exported from E. coli by efflux pumps, predominantly AcrAB-TolC [21]. Before the identification of AcrAB-TolC, Thanassi et al. [37] postulated the existence of and computed kinetic parameters for this pump with tetracycline as a substrate. Later, these kinetic parameters were measured experimentally for ampicillin [21]. The parameters used in our model have been compiled in Table SM6.
We assumed Michaelis-Menten kinetics with potential cooperativity, so the efflux rate v is: In the literature, v max values (nmol/mg dry mass/sec) are reported instead of k cat (1/sec) for these reactions, so we computed k cat as follows, where[E] is the average pump concentration, m is the average dry mass, andv p is the average periplasmic volume (all drawn from a singlecell simulation with an initial seed of 0):

Beta-lactam hydrolysis
The E. coli reference genome used by our models (https://www.ncbi.nlm.nih.gov/n uccore/U00096.3) contains the endogenous beta-lactamase gene ampC that exhibits low and non-inducible expression [26]. We modeled the hydrolysis of ampicillin by AmpC as a Michaelis-Menten reaction (Eq. 4) where the hydrolysis product is inert.

Implementation
Diffusion reactions take place too quickly for the Vivarium engine to accurately integrate using the process time step of two seconds. Therefore, the AntibioticTransportOdeint process used the solve ivp numerical integration function from SciPy [32] to simulate the ordinary differential equations (ODEs) described above over the course of each two-second time step.

Tetracycline-induced changes to gene expression
We modeled tetracycline-induced gene regulatory changes using the TFBinding process that was already in the model [23].

Model description
In our model, tetracycline-induced gene regulation began with the inactivation of MarR by tetracycline, which we modeled as a reversible reaction at chemical equilibrium. In the model, as the fraction of inactive MarR increased, so did the occupation of MarA on the promoters of its downstream regulatory targets. The effect of MarA binding to promoters was tuned to yield comparable mRNA fold changes as those measured for E. coli cells exposed to 1.5 mg/L of tetracycline (Fig. S6B) [40]. Notably, we chose to silence MarA activity in the complete absence of MarR, preserving baseline behavior at the cost of rare (< 5%) delays in tetracycline-induced gene regulation.
In addition to direct regulation by MarA, the ompF gene, which encodes the primary porin for tetracycline influx, is also subject to post-transcriptional regulation that is triggered by tetracycline. Specifically, MarA upregulates expression of micF, a small non-coding RNA that can form duplexes with ompF mRNA and thereby decrease ompF translation [9]. In the absence of experimental data, we assumed that all synthesized micF transcripts immediately form irreversible duplexes with free-floating ompF transcripts and that the resulting duplexes have the same half-life as standalone ompF transcripts. Additionally, we assumed that MarA increased production of micF RNA just enough to sequester nearly all ompF transcripts when cells were exposed to 1.5 mg/L of tetracycline.

Tetracycline binding to ribosomes 4.4.1 Inputs and outputs
Inputs: Outputs: • Free tetracycline concentration

Model description
Tetracycline is a bacteriostatic antibiotic that inhibits protein synthesis by binding to a highly conserved pocket in the 30S ribosomal subunit near the A site, interfering with accommodation of aminoacylated tRNA during polypeptide elongation [15; 18]. Thus, we introduced the TetracyclineRibosomeEquilibrium process to model ribosomal binding as a competition between tetracycline and aminoacylated tRNAs. Since the binding constants for tRNAs and tetracycline are both much smaller than the average concentration of active ribosomes (10 −6 M < 10 −5 M ), we assumed that all ribosomal A sites are bound by one or the other at all times. See Table SM7 for the exact binding constants.
Specifically, given a tRNA-ribosome binding constant of K tRN A , a tetracycline-ribosome binding constant of K T c , c tRN A aminoacylated tRNA molecules, c T c tetracycline molecules, c R,T c tetracycline-bound ribosomes, and c R,tRN A tRNA-bound ribosomes, we know that at equilibrium the following holds: Thus, the fraction f of ribosomal A sites bound by tetracycline at equilibrium is equal to: At the time of writing, vivarium-ecoli lacked a proper mechanistic model for tRNA charging.
Since the average cell in our model reported a total tRNA count about three times the literature consensus for charged tRNAs [12; 17], the simulated tRNA count was always multiplied by a factor of 0.33 before further calculations. Additionally, we assumed that each active ribosome had two bound tRNAs at all times (one in the E site and one in the P site), further decreasing the pool of charged tRNAs competing for A site binding.
Since we only found parameters that describe the equilibrium binding of tetracycline and tRNAs to ribosomes and not their binding kinetics, we used a root finder (root scalar from Scipy [32]) to solve for the equilibrium concentrations of free tetracycline, tetracycline-bound ribosomes, and free ribosomes. The algorithm used to compute this equilibrium is described in Algorithm 1.
Algorithm 1: Tetracycline-ribosome equilibrium Input : c tRN A,tot total tRNA count Input : c T c free cytoplasmic tetracycline count Input : c 70S active ribosome count Input : c 30S free small subunit count Input : c 30S-T c tetracycline-bound small subunit count 1. Estimate the count of free, charged tRNAs from the total count. c tRN A = c tRN A,tot · 0.33 − 2 · c 70S 2. Use root finder to calculate ∆ tetracycline molecules to bind to (or unbind from) 30S ribosomes such that the following equation holds: 3. Distribute the calculated count of ribosomes bound to tetracycline proportionally among 30S subunits and 70S active ribosomes.
Decrease c 70S by c 70S-T c,target and increase c 30S-T c by c 70S-T c,target (inactivated 70S active ribosomes assumed to dissociate and form 30S subunits bound to tetracycline).

PBP binding and inhibition 4.5.1 Inputs and outputs
Inputs: • Periplasmic ampicillin concentration • Count of PBP1A, PBP1B (γ isoform) in the cell • Count of newly produced murein Outputs: • Fraction of PBP1A, PBP1B not bound by ampicillin • Allocation of newly produced murein to pools that are either usable (crosslinked) and unusable (uncrosslinked) for incorporation into the cell wall

Model description
We modeled the ampicillin-mediated inhibition of PBP transpeptidase activity using a Hill equation with no cooperativity. This model was implemented in the PBPBinding process.
We calculated the proportions of unbound, active PBP1A and PBP1B as follows: where θ P BP 1A , θ P BP 1B are the proportion of unbound PBP1A and PBP1B respectively, [Amp] is the concentration of ampicillin in the periplasm, and K A,P BP 1A and K A,P BP 1B are the ampicillin concentrations necessary to inhibit transpeptidation activity by 50% for PBP1A and PBP1B, respectively.
K A,P BP 1B was parameterized directly from literature [6], whereas K A,P BP 1A (specifically for inhibition of transpeptidation) has not been measured to our knowledge. However, ampicillin binding affinities for PBP1A and PBP1B have been measured as being within one order of magnitude [7; 19]. Assuming that the ability of ampicillin to inhibit transpeptidation is proportional to its binding affinity, K A,P BP 1A should also be within one order of magnitude of K A,P BP 1B . Using this plausible range as a guideline, we estimated a K A,P BP 1A for ampicillin inhibition of PBP1A transpeptidase activity that resulted in cell death as expected in singlecell simulations exposed to the 2 mg/L ampicillin (MIC).
Once the proportions of unbound PBP1A and PBP1B were determined, all the nascent murein produced by Metabolism in the most recent time step was divided into two pools, one for crosslinked murein and the other for murein that was not crosslinked due to inhibition of PBP-catalyzed transpeptidation by ampicillin. The count of crosslinked murein was reported to the cell wall process for incorporation while uncrosslinked murein was assumed to be permanently unusable for cell wall synthesis. If a cell had c P BP 1A PBP1A and c P BP 1B PBP1B molecules, the count of crosslinked murein c c and uncrosslinked murein c u was partitioned from the total new murein count c m as follows:

Parameter
Value Description K A,P BP 1A 0.7 µM Ampicillin concentration at which PBP1A transpeptidation activity is halved. K A,P BP 1B 1.27 µM Ampicillin concentration at which PBP1B transpeptidation activity is halved [6]. Outputs: • Updated cell wall lattice • Extension factor above resting length • Cracking

Model description
Cell wall growth was modeled as a function of cell shape, available murein, and active PBPs in the CellWall process. Cell wall state was represented as a 2D lattice on the surface of a cylindrical shell, in which lattice positions represented the average surface area spanned by a peptidoglycan unit. Lattice positions filled with ones represented cross-linked peptidoglycan, while lattice positions filled with zeroes corresponded to gaps where peptidoglycan was not crosslinked into the cell wall.
Because the E. coli sacculus is elastic primarily in the long direction of the cell, we permitted lengthwise stretching of the lattice up to the experimentally determined maximum surface increase, E max . At each 10 second time step, CellWall first attempted to relax the lattice if excess murein and PBPs were available. If current cell dimensions provided by the CellShape process indicated that the cell had grown, the difference between the length of the lattice and the cell length was first expressed in terms of the number of new columns to be added.
After distributing the murein available to be incorporated (calculated by the PBPBinding model) uniformly at random among the columns, the content of each new column was determined by sampling several peptidoglycan strands whose lengths followed a geometric distribution with parameter p strand estimated from literature data. These strands (stretches of ones) were placed end-to-end with single-position gaps (0) in between to populate each newly generated column. Then, N sites = min(# active PBPs, # new columns) insertion sites were chosen uniformly at random along the length of the lattice. New columns were then distributed uniformly at random among the N sites insertion sites such that ≥ 1 column was inserted per site.
The CellWall process then inserted these newly sampled columns at their designated insertion sites, resulting in a "proposed" next cell wall state. The size of the largest hole in the proposed lattice was compared with a critical hole area for lysis from literature (A critical = πr 2 critical ). If the proposed lattice crossed this threshold, the CellWall process first checked to see whether the cell could be saved from lysis by stretching the cell wall further to cover the new cell length. When attempting to stretch, the necessary extension factor was first compared with E max to determine whether stretching was possible. If so, the extension factor E was increased such that the lattice covered the new size of the cell without incorporation of new murein. This alternative proposed lattice was then evaluated again for cracking, since the stretching itself increased the physical size of existing lattice defects. A summary of the steps in cell wall growth is outlined in Algorithm 2.
Once the cell wall had cracked, this was read by the LysisInitiation process. This process then sampled a waiting time (on the order of 3.2 minutes) from a distribution fitted to literature data, representing the time that the cell spends with the inner membrane bulging out through the cell wall. Upon reaching the end of this delay period, Lysis was triggered, resulting in the removal of the cell from the simulation and the spilling of internal ampicillin and beta-lactamase into the environment. These environmental changes were enacted by the ReactionDiffusion process (see 3.3). If a cell survived until division without lysing, the two resulting daughter cells each inherited one-half of the lattice from the mother cell.

Murein adjustment
In the process of creating a physical representation of the cell wall, we noted that the original release of the model produced over two times more murein than necessary to completely cover the surface area of an average cell. As such, we iteratively estimated a corrected homeostatic objective for murein that resulted in approximately no leftover murein at the start and end of a representative cell cycle (seed 0). This new objective was about 2.27 times smaller than the original value and was derived from the assumption that cells do not produce significantly more murein than necessary to maintain cell wall integrity.

Parameter
Value Description p strand 7.68 × 10 −2 Strand extension probability. Fitted assuming a geometric distribution using data from [28; 41]. r critical 20 nm Critical defect radius to initiate lysis [8]. ℓ y 1.03 nm Length of one peptidoglycan unit (along the direction of the strand) [42]. ℓ x 1.4 nm Width of a glycan strand [39]. d x 0.6 nm Typical resting distance between glycan strands. Estimated such that the initial lattice uses all initial murein. Together with ℓ x , consistent with literature distances between centers of adjacent strands [4; 39]. E max 3 Maximum permissible expansion of the cell sacculus [20]. t lysis 192.8 s Mean time to lysis after bulging begins [43]. Table SM9: Parameters for the cell wall process.
Algorithm 2: Cell wall growth Input : m crosslinked murein count Input : c P BP 1A , c P BP 1B PBP1A/B counts Input : θ P BP 1A , θ P BP 1B fraction of PBP1A/B that is active Input : l cell length Input : e cell wall extension factor Input : L cell wall lattice with y rows and x columns 1. Calculate new number of columns in cell lattice.
x new = l e · (l x + d x ) 2. Shrink extension factor if there is excess murein.
if ⌊m/y⌋ > x new − x then x new = x + ⌊m/y⌋ , 1 3. Calculate count of murein allocated to each new column. c syn = min(θ P BP 1A · c P BP 1A + θ P BP 1B · c P BP 1B , y new − y) 5. Distribute newly generated columns among c syn positions in the initial lattice. 6. Record hole sizes, try to stretch if hole too large, crack if still too large.
if max hole size > πr 2 critical then e new = l x(l x + d x ) if e new < E max then if max hole size > πr 2 critical then Cell wall cracked. else Cell wall intact.
else Cell wall cracked.