Minimal biophysical model of combined antibiotic action

Phenomenological relations such as Ohm’s or Fourier’s law have a venerable history in physics but are still scarce in biology. This situation restrains predictive theory. Here, we build on bacterial “growth laws,” which capture physiological feedback between translation and cell growth, to construct a minimal biophysical model for the combined action of ribosome-targeting antibiotics. Our model predicts drug interactions like antagonism or synergy solely from responses to individual drugs. We provide analytical results for limiting cases, which agree well with numerical results. We systematically refine the model by including direct physical interactions of different antibiotics on the ribosome. In a limiting case, our model provides a mechanistic underpinning for recent predictions of higher-order interactions that were derived using entropy maximization. We further refine the model to include the effects of antibiotics that mimic starvation and the presence of resistance genes. We describe the impact of a starvation-mimicking antibiotic on drug interactions analytically and verify it experimentally. Our extended model suggests a change in the type of drug interaction that depends on the strength of resistance, which challenges established rescaling paradigms. We experimentally show that the presence of unregulated resistance genes can lead to altered drug interaction, which agrees with the prediction of the model. While minimal, the model is readily adaptable and opens the door to predicting interactions of second and higher-order in a broad range of biological systems.


Parameter reduction and the bifurcation point
In the following we derive expressions that constrain the parameter space in which the system is bistable. We start by rewriting the steady-state solution with new variables, which will facilitate reduction of parameters. We recast Eq (4) so as to express a ex /2IC * 50 as function of growth rate; this is possible since a ex /2IC * 50 occurs linearly in the equation. At this stage we introduce new variables and rewrite parameters in a dimensionless form: (S1.1a) λ = y × λ 0 , (S1.1b) where IC 50 is concentration needed to halve the growth rate and λ 0 is a drug-free growth rate. These definitions lead to: IC 50 IC * 50 cα = 1 4 α 2 + y 2 − y 3 . (S1.2) To remove the ratio IC 50 /IC * 50 from the equation, we use the expression from Ref. [1] that states IC 50 IC * 50 = 1 2 This relation is obtained by plugging y = 1/2 and c = 1 into Eq (S1.2). Finally, we recast Eq (S1.2) into This equation expresses the relative concentration c as a function of relative growth rate y in which α is the sole parameter. Fitting an implicit function can be challenging; we can estimate α from equating the implicitly calculated derivative of Eq (S1.4) to the derivative of the Hill function with steepness parameter n, when both derivatives are evaluated at c = 1. An estimate for the response parameter is then α ≈ 1/(n − 1).

Bifurcation point
To determine the critical value of the parameter α below which bistable concentration region exists, we consider the stationary points of Eq (S1.4). The derivative reads: (S1.5) Considering the equation 8y 3 − 4y 2 + α 2 = 0, (S1. 6) we can determine values of y where the derivative of Eq (S1.5) is 0. Since the equation is of third order, we expect either three real roots or one real together with a conjugated pair of complex roots. The solutions [found using Mathematica (Wolfram Research, 11.3) function Solve] are: We note that whether all roots are real is determined by the sign of 27α 2 − 4; this leads to a critical value of α crit = 2/(3 √ 3) ≈ 0.385-the system has a bistable region when α < α crit . To obtain the span of concentrations at which two bistable growth rates exist (Fig 1D), we plug Eqs (S1.7a) and (S1.7c) into Eq (S1.4).

Linearized case of constitutively expressed resistance gene
As noted in the main text, for high values of K rem the Michaelis-Menten equation linearizes. In this case we can reuse the expression for growth-dependent efflux from Ref. [1]. Using the parameter reduction from above, the generic solution is where Ω = V max /(p out K rem ). This expression can be rewritten into reduced form as Here, we used the relation that follows directly from Eq (S1.8); IC 50 /IC * 50 = (1/2α) × 1 + α 2 (1 + Ω/2) . By solving dc/dy = 0 we obtain a condition −64α 2 + 432α 4 + 48α 4 Ω − 12α 6 Ω 2 + α 8 Ω 3 ≤ 0 for bistability to exist; heuristically, we note that the leading term in Ω is of the third power and positive and will for high values of Ω make the expression necessarily positive and will thus make the system monostable. This illustrates the importance of non-linearity in the expression describing antibiotic removal [Eq (19)]; the presence of positive feedback alone is not sufficient for the occurrence of bistability. In the limit Ω → ∞, Eq (S1.9) becomes simply y = 1 − c/2 for c ≤ 2 and y = 0 otherwise.

Calculation of dose-response surface for additive interaction and Loewe interaction score
An additive interaction is characterized by linear isoboles. In this section we briefly sketch the derivation. Let c A and c B be relative concentrations (measured in IC 50 of respective antibiotics) of antibiotics A and B, respectively. Individual dose-response curves are given by f A and f B . To construct the additive surface, the growth rate in any point (c A , c B ) has to be calculated solely from known responses to individual drugs. As isoboles are linear, there is an isobole that passes through this point but terminates in some unique (c A,0 , 0) and (0, c B,0 ), from which it follows that (9) is performed over an area where growth rate is above a chosen threshold. Throughout the article this threshold is set to y = 0.2. Positive or negative values of LI suggest antagonistic or synergistic interaction, respectively.

Numerical solutions
We evaluated the steady-state solution of the system Eqs (6) by forward time integration of the differential equations. We used λ 0 = 2 h −1 , which is the drug-free growth rate in rich lysogeny broth medium at 37 • C [2]. We rescaled the time by defining To mimic the situation in which exponentially growing unperturbed bacterial culture is exposed to antibiotic, we set the initial condition to y = 1; the rest of the species were set to 0. We integrated the rescaled differential equations forward in time until τ = 9 × 10 3 by using Mathematica function NDSolve. In computation, K D and k on values were set for both antibiotics to 0.1 µM and 100 µM −1 h −1 , respectively. The results are largely invariant of the choice of these numeric constants as long k on κ t and K D is roughly between 0.01−1 µM. 4 Analytical solution in the limit of strong inhibition and reversible binding The system of ODEs [Eqs (6)] can be linearized for near-zero growth rates (i.e., λ λ * 0 , λ 0 ) which is the case when the external concentrations of antibiotics are high enough (a ex,i IC 50,i ). In the latter case one can postulate two constraints: (i) Ribosome synthesis is up-regulated to the theoretical maximum, i.e., r tot = r max and (ii) internal concentration of the i-th antibiotic is a i ≈ a ex,i p in,i /p out,i . These constraints eliminate the dynamical equations for internal concentrations of antibiotics [Eq (6a)]. The fully specified and expanded system of equations readṡ The system is linear and thus the steady-state solution exists in a closed and unique form. We initially derive the stationary solution for independent binding, i.e., δ σ,i = 1: Each product term can be rewritten as a function of a i /K D,i , which can be further transformed by noting that K D,i p outi = (λ 0 α i ) 2 / (4κ t ) and a ex,i p in,i = c i × ∆rλ 0 α 2 i + 1 /4. Together, these expressions allow rewriting Eq (S4.11a) using variables c i = c i × (α 2 i + 1)λ max /(α 2 i λ 0 ), which finally leads to the Eq (10). To obtain Eq (11) the general solution of Eqs (S4.10) is evaluated; the latter is further simplified by assuming that k off,A ≈ k off,B , which approximately holds for reversibly binding antibiotics.

Effect of dose-response curve concavity on the shape of isoboles
For intermediate values of α ∼ 1 and increasing δ, we observed that the isoboles of the dose-response surface at lower drug concentrations indicate strong antagonism, whereas at higher concentrations, they indicate synergism (Fig S1A). Is it possible to determine a concentration value above which increasing δ For α < 2 black dashed line denotes the boundary below which the increasing cooperativity δ increases antagonism (shades of yellow); above it synergistic character is enhanced (blue). Two isoboles are showcased for each example: purple and black lines correspond to δ = 1 and δ = 10 4 , respectively. Note, that below c inf isoboles in the case of strong cooperativity become steeper than the independent ones; above c inf isoboles become shallower, which is indicative of synergy.
will cease to increase antagonism but rather increase synergism? Intuitively, the interaction-characteristic shape of isoboles is determined by the sign of the mean curvature, the latter being defined as where k 1 and k 2 are principal curvatures. Full dependency of K a is not accessible as an analytical expression for y(c A , c B ) is not known; however in the limit δ → ∞ the antagonistic isoboles become nearly perpendicular to the axes, thus rendering all derivatives along perpendicular directions equal to zero. Hence, the expression for K a is simplified into: As the denominator is always positive the sign of K a depends only on the sign of the second derivative along the axis. To find this point, we need to determine the inflection point of the relative growth rate y along either of the axes. Here, we need to solve d 2 y/d 2 c = 0 for α; since y is given only implicitly in Eq (5), the second implicit derivative is calculated. This leads to inflection point at relative growth rate y inf = 3 α 2 /4 for α > α crit , at which second derivative vanishes. Using Eq (5) the inflection concentration is evaluated as: (S5.14) From the expression for y inf we observe that we have a relevant solution only for α < 2; this leads to the first conclusion-the intermediate regime of surfaces having both antagonistic and synergistic contours exists only up to α = 2. This also implies that for every response parameter α, given a high enough concentration of an antibiotic, the isoboles will become synergistic. However, at low growth rates at high concentrations the approximation [Eq (13)] becomes relevant (Fig S1B). Numerically computed doseresponse surfaces for δ = 1 and δ = 10 4 illustrate that indeed above the c inf character of the isoboles is different and that the transition happens due to vanishing second derivative.

Experimental considerations
We verified the results on antibiotic induced starvation and constitutive expression of resistance genes by measuring the growth rate in two-dimensional gradients as described in Ref. [2]. Briefly, concentration gradients were constructed in the wells of microtiter plates; multiple such plates were cycled through the plate reader, which measured the luminescence as the proxy of the cell density. From the time traces of luminescence the growth rate of the bacterial culture in the individual well was determined by identifying an exponential region and fitting the line to the log-transformed luminescence values-precise methods and experimental setup description is in Ref. [2]. Bacterial growth medium used in all experiments was lysogeny broth (LB). Antibiotics solutions were prepared from powder stocks and dissolved in ethanol (CHL, MUP), water (STR) or directly in LB [TET and CHL (for CERG experiment)] and filter-sterilized. Growth rates were normalized with respect to the antibiotic-free growth rate λ 0 to obtain y = λ/λ 0 ; concentrations were normalized with respect to IC 50 , which was determined from fitting a Hill function 1/ [1 + (a ex,i /IC 50,i ) n ]. To construct the double-resistant strain we followed the methods described in Refs. [2] and [3] with appropriate modifications. We have cloned cat gene into a plasmid with pKD13 background (as per Ref. [2]), in which resistance gene was driven by a synthetic promoter P LlacO−1 [4]. Promoter was unregulated (constitutive) as the background strain HG105 (MG1655 ∆lacIZYA) [5] is devoid of the entire lac operon, including the lac-repressor. We amplified the tetA and cat from the strain MS004A [6] and plasmid pZA32 [4], respectively. We were unable to clone tetA into a plasmid; we therefore assembled the kanamycin cassette, P LlacO−1 promoter, and tetA gene with a rrnB terminator in vitro using HiFi Assembly Mix (NEB), PCR-amplified the fragment, and integrated it into the intS locus. CAT-coding gene was integrated into galK locus. Both CERGs were integrated into the chromosome by lambda-red recombineering, selected on kanamycin, which was followed by the resolution of the kanamycin cassette by FLP-resolvase [7,3,8,2,9]. We verified modified chromosomal segments with PCR and Sanger sequencing.

Combinations of translation inhibitors with drugs that alter growth law parameters
We elaborated in the main text on a specific example of a starvation mimicking antibiotic, which invokes a response that is similar to shifting bacteria to poorer media (i.e., change in λ 0 ). As growth laws capture this change, a response could be analytically predicted and was experimentally tested. Below we show two more hypothetical examples of antibiotic combinations in which translation inhibitor is combined with a drug that alters growth law parameters. We consider the case in which this hypothetical drug alters either r min alone or in concert with r max ; in either case a shift in growth laws happens, which affects both dynamic range ∆r as well as the apparent response parameter α. Here, we note that there is some experimental data illustrating that such shift could be achieved by antibiotics inhibiting transcription, Purple arrow denotes the effect of growth law-varying drug, which alter the "origin" of the translation inhibition growth law line (dashed gray lines). Slope of dashed lines is altered when r min is varied; when r max is varied simultaneously with r min the slope is invariant of the "origin." Purple squares is data describing the ribosome abundance in response to rifampicin from Ref. [10]; we converted reported RNA/protein measurements to ribosome concentrations by fixing the drug-free data point to the ribosome concentration r u ≈ 32.0 µM as determined by the Eq (1) at λ 0 = 0.78 h −1 . Note, that the data points lie approximately on the horizontal dashed line. (B) Examples of dose-response surfaces for different translation inhibitors (concentration c; different response parameters α) and different dose-response curves for ∆r-and r min -varying antibiotic (concentration x; Hill-function with different steepness parameter n). For easier comparison we overlaid the dose-response surfaces; the purple contours correspond to ∆r-varying case.
e.g., rifampicin (Fig S2A and supplementary information of Ref. [10]). However, the combined effect of transcription and translation inhibition of translation machinery has not be experimentally assessed.
If we consider an antibiotic that affects the r min , yet it does not affect the κ t (Fig S2A), then the following modification to the growth law arises: Here, we denoted by g min (x) a dose-response function of a r min -varying antibiotic and x its concentration; r min denotes an apparent new r min as a function of g min (x). Thus, the reduced dynamic range reads ∆r = r max − r min = ∆r − [1 − g min (x)] λ 0 /κ t = ∆r {1 − [1 − g min (x)] r 0 /∆r}, where we defined r 0 = λ 0 /κ t . Taken together, these relations rescale the response parameter α = α F /g min (x) (where α F is the response parameter for x = 0), while the apparent concentration becomes c = c/ψ min , where ψ min = IC 50 IC 50,F = (S7.16) = α 2 F + g 2 min (x) g min (x) (α 2 F + 1) If we consider that g min (x) = 1/ (1 + x n ), where x is the concentration of r min -varying antibiotic measured in units of IC 50 , we can calculate the whole dose-response surface as y(c, x) = g min (x) × f (α F /g min (x), c/ψ min ) , (S7.17) where f is a translation inhibitor dose-response curve. Examples are shown in Fig S2B. Model variation discussed above can be expanded in a model in which both r min and r max are varied simultaneously (Fig S2A). Here, we assume that r min increases and r max decreases in concert such that they have the same value at g ∆r = 0 (Fig S2A). It follows ∆r = r max − r min = = r max g ∆r (x) + [1 − g ∆r (x)] (r 0 + r min ) The response parameter α is rescaled as above; the ratio ψ = IC 50 /IC 50,F takes a slightly different form: With this at hand, dose-response surfaces can be evaluated (Fig S2B). Here we note that the both models give the same result if r 0 = ∆r, which requires λ 0 = ∆rκ t . This is intuitive since the ribosome concentration is already at its maximum and only r min can be varied.