Individual-based model highlights the importance of trade-offs for virus-host population dynamics and long-term co-existence.

Viruses play diverse and important roles in ecosystems. In recent years, trade-offs between host and virus traits have gained increasing attention in viral ecology and evolution. However, microbial organism traits, and viral population parameters in particular, are challenging to monitor. Mathematical and individual-based models are useful tools for predicting virus-host dynamics. We have developed an individual-based evolutionary model to study ecological interactions and evolution between bacteria and viruses, with emphasis on the impacts of trade-offs between competitive and defensive host traits on bacteria-phage population dynamics and trait diversification. Host dynamics are validated with lab results for different initial virus to host ratios (VHR). We show that trade-off based, as opposed to random bacteria-virus interactions, result in biologically plausible evolutionary outcomes, thus highlighting the importance of trade-offs in shaping biodiversity. The effects of nutrient concentration and other environmental and organismal parameters on the virus-host dynamics are also investigated. Despite its simplicity, our model serves as a powerful tool to study bacteria-phage interactions and mechanisms for evolutionary diversification under various environmental conditions.

Introduction

Viruses execute a wide range of functions in the biosphere, influencing biogeochemical cycles, affecting efficiencies of transport of energy and matter through food-webs and driving processes of evolutionary diversification [1-4]. These ecosystem-related functions may seem disentangled from the effects that pathogenic viruses can have on human health and society as demonstrated by the ongoing COVID-19 pandemic [5]. Yet, the underlying mechanisms driving the dynamics between viruses and their hosts are in principle the same.

Viral ecology has been growing as a field of research in the last decades, with sequencing techniques revealing enormous biodiversity in viral genomes [6]. The field, however, remains challenging, since advanced and at times intricate laboratory experiments are required to characterize viral traits as well as interactions with their hosts. Besides, much of the environmental viral metagenome remains unmapped or still undiscovered [7]. We conjecture that virology will grow as a field if we manage to focus on principles that unify different disciplines, from viral ecology to infectious disease research. To better understand the fundamental mechanisms that drive viral dynamics, conceptual models should be used as tools to identify principles that explain biodiversity and functioning in viral systems, be they ecological feedback mechanisms or emerging evolutionary dynamics.

Trait-based approaches have proved useful in identifying unifying and universally applicable mechanisms in ecology and were first established in terrestrial ecology [8, 9]. They have also been successfully applied to marine ecology [10], and in particular to microbial ecology [11-13]. A strength of trait-based approaches is that processes and interactions between organisms and their environment are described on a functional level, independent of particular taxa at hand. Besides, it brings trade-offs between organismal traits to the center of attention, which arise from chemical and physical constraints that are universal (ie, system-independent). Such trade-offs are crucial to understanding life [14-18]. In the marine microbial ecosystem, various processes and structures have been linked through fundamental trade-offs [19]. Trait- and trade-off based perspectives thus also promise to unify viral ecology [20].

In the present context, we use a simple evolutionary individual-based model to test the hypotheses that trade-offs between competitive and defensive traits are key to understanding virus-host population dynamics and evolutionary change. We define competitive traits as organism-specific traits to acquire limiting resources (in our case, nutrient affinities of hosts) and defensive traits as traits modifying the efficiency of a viral infection (here, the viral adsorption coefficient and the inverse of the host range of viruses). Following [21], we model each of these traits as genes being embedded in a phenotypic trait function, and we incorporate trade-offs between them by means of interaction functions that determine the likelihood of infection based on gene similarity of the host and virus (“compatibility function”) and virus-intrinsic infection efficiency (“virulence function”). The model is highly idealized, with emphasis on infection mechanisms that are important and applicable to both environmental as well as pandemic settings.

In the next section, we describe our individual-based model along with experimental methods used to characterize virus-host infection dynamics. This is followed by findings from our model, including a validation of the model dynamics with experimental results. We conclude the article by discussing the relevance and generality of our results for increased insights into virus-host interactions and evolutionary dynamics.

Materials and methods

Individual-based model

Individual-based models (IBMs) are in-silico models that describe the behavior of autonomous individuals (organisms). These models are widely used, not only in ecology [22] but also in other disciplines dealing with complex systems made up of autonomous entities [23]. In this section, we give an informal description of an evolutionary IBM to study the interaction patterns between bacteria and virus. For a formal description of the algorithm used in the IBM, we refer the reader to the pseudocode for the IBM in the Supporting Information section. An implementation of the IBM in the Python programming language is available at a public Github repository [24].

The state variables and parameters (environmental and organismal) of the IBM are listed in Table 1, each with a symbol, description, typical value [21], and units.

Table 1

State variables and parameters of the IBM.

Symbol	Description	Value	Units
State variables
N h	host abundance	variable	individuals
N v	virus abundance	variable	individuals
Environmental parameters
T	simulation duration	3.6 × 101	h
t	time step duration	1	h
V	chemostat volume	10−6	L
P	total phosphorus concentration	7.1 × 101	μmol-P L−1
P d	dissolved phosphorus concentration	variable	μmol-P L−1
ω	chemostat dilution rate	2 × 10−1	h−1
Nh0	initial value for Nh	4.6 × 101	individuals
Nv0	initial value for Nv	8.1 × 102	individuals
VHR	virus to host ratio (the ratio Nv0/Nh0)	1.76 × 101	-
Organismal parameters
P h	maximum phosphorus concentration in a host	8.3 × 10−8	μmol-P
α	nutrient affinity of a host	1.6 × 10−7	L h−1
μ h	maximum growth rate of a host	7.38 × 10−1	h−1
g α	nutrient affinity gene of a host	variable ∈ [0, 1]	-
gα0	initial value for gα	10−1	-
m	mass of a host	variable ∈[12,1]	-
π h	probability of mutation of host genotype	6 × 10−2	-
σ h	standard deviation of host genotype mutations	5 × 10−2	-
δ h	mortality rate of a host	1.4 × 10−2	h−1
ϵ h	metabolic loss rate of a host	1.4 × 10−2	h−1
β	adsorption coefficient of a virus	6.2 × 10−11	L h−1
g ν	memory gene of a virus	variable ∈ [0, 1]	-
gν0	initial value of gν	10−1	-
g β	adsorption coefficient gene of a virus	variable ∈ [0, 1]	-
gβ0	initial value of gβ	10−1	-
π v	probability of mutation of virus genotype	6 × 10−3	-
σ v	standard deviation of virus genotype mutations	5 × 10−3	-
κ	number of viruses produced per infection of a host	101	-
δ v	decay rate of a virus	1.4 × 10−2	h−1

Symbol, description, typical value [21], and units.

For simplicity, a single elemental resource is used in the model budget. Specifically, phosphorous is used as the model currency [25]. For ease of interpretation of the results, the total phosphorus concentration P is expressed in terms of the number of host individuals in the simulated volume, ie, P ∼ P/(P ⋅ V). Similarly, the nutrient affinity α of a host and the adsorption coefficient β of a virus are normalized in terms of the simulated volume, ie, α ∼ α/V and β ∼ β/V.

At the start of the simulation, we consider a nutrient medium of volume V. The medium is inoculated with a host population , in which each host h has genotype and mass m picked uniformly at random from the interval . The medium is also inoculated with a virus population , in which each virus v has genotype . Note that H and V are sets of individuals, which grow (or shrink) as dynamics unfold. The initial amount of dissolved phosphorus is calculated in units of host individuals as .

At each time step t ∈ [1, T] of the simulation, we update the current concentration P of dissolved phosphorus as , which takes into account the inflow and outflow of phosphorus due to washout. We then carry out a round of host dynamics followed by a round of virus-host interaction dynamics as described below. Each time step t is fixed to 1 hour. With time-steps lasting 1 hour, rates in units h−1 can be translated into probabilties per time step.

During the host dynamics round, we consider each host h ∈ H in sequence. The host can be washed out of the medium with probability ω. If it remains in the medium, it can die with probability δ, in which case its biomass m is returned to the medium, ie, P is incremented to P + m. If the host does not die, it experiences a mass loss due to metabolism given by ϵm, which is subtracted from m and added to P. The host h grows during the time step t, and the resulting gain in mass is calculated as [26], which is subtracted from P and added to m. If m exceeds unity, then the host divides into two daughter cells, each with half the mass (ie, m/2) and the same genotype (ie, {g}) as the host. With probability π, the genotype of the daughter cells mutates to , where is sampled from a Gaussian distribution with mean g and standard deviation σ; and with probability 1 − π, the genotype of the daughter cells remains the same as that of the parent cell. The daughter cells are added to the host population H and the parent cell is removed from it.

During the virus-host interaction round, we consider each virus-host pair (v, h), where v ∈ V and h ∈ H. The virus v can be washed out of the medium with probability ω. If it remains in the medium, it can decay with probability δ. If the virus does not decay, then whether or not it infects the host h depends on its compatibility with the host and its virulence. The former is given by the compatibility function and the latter by the virulence function . The values of both functions (also called interaction functions) are from the unit interval [0, 1], and are interpreted as probabilities. To test the effects of trade-off based vs random compatibility between hosts and viruses, we consider two compatibility functions: trade-off based compatibility and random compatibility , which samples a number uniformly at random from the unit interval. The function (Fig 1) captures the trade-off between the host’s nutrient affinity and the virus’ adsorption coefficient, such that the compatibility is highest when the corresponding gene values g and g are similar. The virus-host compatibility is also enhanced when the virus’ host range (expressed by its memory gene g) is high.

Fig 1

Trade-off based compatibility function .

The left panel shows how varies with |g − g| for fixed values of g; is high when g and g have similar values and low otherwise. The right panel shows how varies with g for fixed values of |g − g|; is high for high values of g and low otherwise.

We also consider two virulence functions: trade-off based virulence and random virulence , which samples a number uniformly at random from the unit interval. The function (Fig 2) captures the trade-off between the virus’ adsorption coefficient g and its host range g, whereby virulence of the infection is high for viruses with high adsorption (ie, high g) and narrow host range (ie, small g). The virus v is compatible with host h with probability . If they are compatible, the virus can infect the host with probability . All viruses are assumed to by lytic, meaning that hosts die upon infection. If the host is infected, then it is removed from the host population H, its biomass is immediately recirculated into the dissolved nutrient pool P, and the virus produces κ copies of itself, each of which has the same genotype (ie, {g, g}) as the parent. With probability π, the genotype of the copies mutates to , where and are sampled from a Gaussian distribution with mean g and g respectively and standard deviation σ; and with probability 1 − π, the genotype of the copies remains the same as that of the parent. The virus v cannot infect any more hosts. It is removed from the virus population V and the copies of the virus are added to the population.

Fig 2

Trade-off based virulence function .

The left panel shows how varies with g for fixed values of g; is high for high values of g and low otherwise. The right panel shows how varies with g for fixed values of g; is high for low values of g and low otherwise.

The behavior of the model is investigated for different parameter values and interaction functions. Sensitivity analyses for virus to host ratios (VHR), limiting resource concentration, and dilution rates are carried out to study the effects of environmental conditions on virus-host dynamics. A range of host and virus mutation probabilities as well as standard deviation for the mutations are considered to see how the arms-race dynamics are influenced by host- and virus-specific cellular constraints. Finally, genotypic and random compatibility/virulence functions are explored in order to analyze the sensitivity of arms-race dynamics to trade-off based virus-host interactions. Additional investigations of the effects of various physiological parameters on virus-host population dynamics are summarized in the supplementary material. In all cases, the effects of different parameters and interaction functions are ascertained from an ensemble average of 100 independent simulations.

Laboratory experiment

To characterize virus-host infection dynamics in a biological system, as well as to validate the general dynamics emerging from our IBM, we performed an infection experiment with the bacteria Escherichia coli E28 (DSM 103246) and a T-4 like virus (DSM 103876, hereby called B28) at five different VHRs (0.01, 0.15, 0.73, 1.90, 3.30). The experiment was conducted in a 96-well flat-bottomed microplate using the 2300 EnSpire Multilabel Plate Reader (PerkinElmer), allowing the entire experiment including multiple replicates to be run simultaneously under one assay. The growth of the host was monitored through automated measurement of optical density at 600 nm (OD600) every 15 mins. All cultures were grown in the minimal medium M9 containing 2 mM MgSO4⋅7H2O, 0.1 mM CaCl2⋅2H2O, 6 mM glucose [27]. All dilutions were also done using M9 medium as diluent.

An outline of the steps involved is as follows: an overnight culture of E. coli was prepared, adjusted to OD600 = ∼0.2 using a Cell Density Meter (Fisher Scientific), and further diluted by factor 1:25. A sterile microplate was filled with the host culture, as well as M9 medium as blank solution, final volume of 200 μL per well. The assay was then run at 37°C, 150 rpm (linear mode, 3 mm diameter) inside the plate reader. The growth of the host culture was monitored and once the host had entered the exponential growth phase, and its OD600 had increased by 0.04 (∼3.5 hr), the assay was paused to retrieve the plate. A sample of the host culture was withdrawn and flash-frozen in 20% glycerol for host enumeration at B28 infection timepoint. B28 lysate (∼108 PFU mL−1) was then added to the host culture to final volumes of 200 μL and concentrations of 0.1, 1, 5, 12.5, and 20% v/v. A sample of the lysate was also flash-frozen in 20% glycerol for phage enumeration. The assay was then resumed under the same culture conditions for another ∼37 hr. Host and phage enumeration was later done using a Calibur flow cytometer (Becton Dickinson) following a standard protocol [28], in order to calculate the exact VHRs. The experimental data (experimental_data.csv) are available as supporting information.

Results

In this section, we present host-virus infection dynamics from our growth experiments in the lab (Fig 3), along with results from our IBM (Figs 4–6). The latter are based on an ensemble average of 100 simulations. First we present the effects of environmental conditions on the outcome of the virus-host interaction dynamics in the first 1.5 days after virus infection (Fig 4), then we show the effects of evolution at the cellular level—specifically, mutation probability and mutation variance—(Fig 5), and finally we show the effects of trade-offs in host-virus interactions, considering a longer time-scale of 30 days (Fig 6). All parameter values except those being tested were held constant as shown in Table 1. Parameter values being tested (Figs 4 and 5) are shown in the figure legends. In all figures, host and virus population dynamics are reported in the first and second row, respectively. Note that in each figure, for ease of comparison, we use the same scale for the y-axes for hosts and viruses.

Fig 3

Infection dynamics of B28 virus and E. coli.

Infection experiments performed as plate reader assays for different virus to host ratios (VHR). Optical density at 600 nm (OD600) was monitored to serve as a proxy for the host abundance. The growth curve of E. coli without virus infection (control) is shown in red. Plotted curves: mean values (n = 8 and ncontrol = 11) shown as solid lines and standard deviaitons shown as shading, with the vertical dashed line denoting the time (t = 3.5 hr) for viral infection.

Fig 4

Model’s sensitivity to environmental parameters.

Population dynamics for hosts (a-c) and viruses (d-f) for varying virus to host ratios (VHR), total phosphorous content (P, μmol-P L−1) and washout rate (ω, h−1). Plotted curves: ensemble averages from 100 runs lasting T = 36 h, with standard deviations shown as shading.

Fig 6

Model’s sensitivity to trade-off-based vs random infection processes.

Host (a-d, i-l) and virus population dynamics (e-h, m-p) shown for closed systems and open systems with washout rate 0.2 h−1 for genotypic compatibility and genotypic virulence (, , upper left), genotypic compatibility and random virulence (, , upper right), random compatibility and genotyipc virulence (, , lower left) and random compatibility and random virulence (, , lower right). Ensemble averages of 100 simulations (lines) with standard deviation (shaded areas) lasting for T = 720 h are shown. Sub-panels show virus dynamics on adjusted y-axis scales.

Fig 5

Model’s sensitivity to evolution in organism traits.

Sensitivity of host (a-d) and virus dynamics (e-h) to mutation probability in host traits (π), virus traits (π), standard deviation in mutation for host traits (σ) and standard deviation in mutation for virus traits (σ). Plotted curves: ensemble averages from 100 runs lasting T = 36 h, with standard deviations shown as shading.

Growth experiments with different VHR over the course of 1.5 days reveal that the host population crashes fastest for highest VHR, thereby also reaching a lower maximum population size at the initial population peak. At the same time, the host population recovers fastest for high VHR, resulting in a pronounced second population peak after one day of incubation (yellow, pink and green curves vs black and blue curves in Fig 3). This has the practical consequence that high VHR allows us to observe a full virus-host infection cycle with subsequent recovery of the host within the timeframe of our experiment.

Our IBM captures these dynamics (Figs 4 and 5). Specifically, analogous to the laboratory experiment, high VHR results in earliest crash but also fastest recovery and emergence of potentially resistant hosts in the infected population within the first 36 hours of our in-silico experiments. This manifests as a second peak in host populations reaching between 800 and 1,100 cells μL−1 roughly 20 hours into the infection cycle (Fig 4a).

Virus numbers increase as they are released from infected hosts and the host population starts collapsing, roughly 10 to 20 hours into the infection cycle. They keep increasing after recovery of the host population, reaching up to 70,000 cells μL−1. Highest VHR at the start of infection cycle consistently results in higher virus population numbers throughout the simulated time frame (Fig 4d).

Nutrient concentration also influences population dynamics; our model indicates that high availability of nutrients leads to a more pronounced boom and bust scenario in the host population, with hosts reaching higher maximum population numbers in their first peak before they crash (up to 16,000 cells μL−1) and growing back to high population numbers (around 1,000 cells μL−1) due to potentially resistant hosts emerging (Fig 4b). Higher limiting nutrient availability is also reflected in highest virus population numbers, which exceeds 80,000 viruses μL−1 after 36 hours (Fig 4e).

Washout also has a strong effect on population dynamics. A crash in the host population with recovery half-way through the simulated time is most pronounced at low washout rates, reaching numbers down to 300 cells μL−1 around 18 hours into the infection cycle. Highest washout rates yield no reduction in host population after viral addition, resulting in steady host population numbers of around 1,200 cells μL−1 from about 12 hours into the infection cycle (Fig 4c). Virus population numbers are highest for low washout rates (Fig 4f). Interestingly, a saturation of the medium with hosts and absence of viruses occurs when the washout value exceeds 0.3 h−1 (Fig 4c and 4f).

Besides the effects of environmental parameters (Fig 4), our model shows clear dependence of the virus-host dynamics on the organism-specific evolution of traits, namely mutation probability π and standard deviation σ of mutations (Fig 5). In particular, increasing the values for the host (π and σ, respectively) leads to a reduced crash and earlier recovery of potentially resistant hosts, with the host population growing back to high densities up to 900 cells μL−1 within the simulated time (Fig 5a and 5c).

For the lowest tested π, hosts struggle to recover within the simulated time, whereas they collapse entirely for the two lowest tested values σ (Fig 5a & 5c). Virus population numbers reach high values of up to 80,000 viruses μL−1 in the simulated time for high π and σ, increasing readily as the host population recovers from the first crash, whereas low π and σ result in stagnant or declining virus population numbers from about 12 hours into infection cycle, never exceeding 40,000 viruses μL−1 (Fig 5e and 5g).

In contrast to host mutation probabilities, changing the mutation probability for the virus π over two orders of magnitude does not show any effect on the host-virus dynamics (Fig 5b and 5f). Similarly, increasing the standard deviation for virus mutations σ does not show a clear effect on virus population dynamics within the simulated time, but the host population recovery after the first crash is dampened with higher σ, reducing the size of the second peak in the host population reaching moderate numbers of 300 cells μL−1 roughly 28 hours into infection cycle and also allowing for a second crash of the host population within the simulated time of 36 hours (Fig 5d and 5h).

Long-term simulations of 30 days (720 h) reveals marked differences of trade-off based vs random interaction function for virus-host compatibility and virulence (Fig 6). Recall that a trade-off based compatibility function () in our simulations implies that infection success is mediated by a trade-off between the host’s nutrient affinity and the virus’ adsorpiton coefficient (Fig 1), whereas a trade-off based virulence function () implies that the infection success is dictated by a trade-off between the virus host range and the virus adsorption coefficient (Fig 2).

In all long-term simulations (Fig 6), the system tends towards an equilibrium after the first fluctuations of host crash and typical re-growth thereof that is also visible in the short-term simulations (Figs 4 and 5). The initial crash in host population is most pronounced in closed systems, and hosts as well as viruses die out completely in closed systems under random compatibility modes (Fig 6i, 6m, 6k and 6o).

In simulations with compatibility trade-offs (Fig 6a–6h), hosts reach a high population size at equilibrium after initial drops both in closed and open systems. As in case of random compatibility mode, the first crash of host population is most pronounced in closed systems (Fig 6a and 6c). Viruses reach very high population size at first, independently of whether its virulence is trade-off based or random, but they experience a steady decline over the course of the simulation. Note that the spread in virus population size over the 100 simulation runs is large in closed systems with trade-off based compatibility compared to all other scenarios. In contrast to closed systems, open systems (washout 0.2 h−1) suppress the emergence of large viral population all together, keeping them at low numbers after an initial spike when simulations are driven by trade-off based virulence (Fig 6f). For open systems with random virulence, virus population size drop right from the start of the simulation (Fig 6h).

In case of random compatibility (Fig 6i–6p), the host population typically collapses to extinction with the consequence that also viruses eventually die out. The only exception is in open systems where virulence is random. In this case, the viruses disappear early on in the simulation and the hosts persist with high population numbers (Fig 6l). With random compatibility and trade-off based virulence in an open system, host population decreases more slowly and undergoes fluctuations while doing so (Fig 6j). Highest viral population numbers before decline are reached in closed systems. As in open systems with trade-off based compatibility but random virulence (Fig 6h), virus population numbers drop from the beginning of the simulation for open systems and random virulence (Fig 6p).

Discussion

We have developed a simple individual-based virus-host interaction model to i) study the effects of environmental and organism-specific parameters on population dynamics and ii) evaluate the importance of trade-offs between competitive and defensive host traits on population dynamics and long-term evolutionary change. The model is highly idealized, focusing on specific trait-based mechanisms. Whereas experimental support for trade-offs between competition and defense in diverse biological systems exists [14, 15, 18], other mechanisms such as environmental disturbance and resource specialization also promote diversity. The model is validated with laboratory data from a short-term infection experiment. Some quantitative discrepancy appears between dynamics in the laboratory experiment and our simulations for low VHR, possibly because all viruses in our model are infectious, whereas a fraction of viruses is non-infections in the laboratory. This implies that higher VHR in the laboratory functionally correpsond with lower VHR in our simulations. Regardless of this, the qualitative match between laboratory data and simulation results for host dynamics provides confidence in basic mechanisms underlying virus—host dynamics. Besides, the model produces plausible results in terms of compatibility trade-offs promoting long-term co-existence of hosts and viruses. However, further validation of the model is needed, in particular for viral population dynamics, which are not as straight forward to measure in laboratory experiments. In the following we discuss details of our findings.

Our model suggests that high VHR results in earliest crash but also fastest recovery of the host population, such that a full infection cycle is complete within the first 36 hours of our simulation experiments. This agrees qualitatively with infection cycles observed in our E. coli-B28 experimental model system (Fig 3) where highest VHR also led to fastest regrowth. Interestingly, similar dynamics have been observed previously in marine algal host-virus systems, namely Emiliania huxleyi and EhV-99B1, Pyramimonas orientalis and PoV-01B, and Phaeocystis pouchetii and PpV-01 [29]. Considering selection pressure to be the highest under strong viral control, it is reasonable to assume that re-growth of host populations after the first crash is given by resistant mutants, which readily establish in high VHR conditions. Virus numbers increase simultaneously with the first crash in host populations as they are set free from infected hosts, and keep increasing after recovery of the host population. This suggests that regrowth of the host population is either not exculsively caused by resistant types, but that it is instead diversifying into resistant and susceptible hosts, or that viruses undergo evolutionary change establishing strains that are able to re-infect the originally resistant hosts. The two alternatives not being mutually exclusive, genomic studies will be needed to decipher the cause.

Analogous to microbial blooms under favorable growth conditions, high availability of limiting nutrient gives a boom-and-bust scenario in our simulations with host population reaching very high peaks before crashing to values below those of more stable host population numbers at lower nutrient concentrations (Fig 4b).

Washout rate is also an important parameter for virus-host population dynamics, especially in shorter time-scales, where washout seems to prevent the accumulation of high viral numbers, keeping viral pressure relatively low in open systems (Fig 6, closed vs open system runs). A possible explanation for low viral pressure in open systems could be a combination of washout of viruses and favorable growth conditions for hosts with inflow of fresh medium.

Evolution of organism-specific traits also affects virus-host population dynamics. Our results suggest that variation in host traits may have stronger effects on virus-host population dynamics than variations in virus traits (ie, different host mutation probabilities and standard deviations give more distinct results compared to different viral mutation probabilities and standard deviations, Fig 5). Noticeably, the host population recovery after the first crash is dampened the most with higher standard deviation in mutation for viruses (σ), reducing the size of the second peak in the host population and also allowing for a second crash of the host population within the simulated time (Fig 5d). In other words, high σ and thus higher flexibility in adjustments of infectious traits as adaptation to acquired host immunity seems to impose a strong selection pressure on hosts, speeding up evolutionary time-scales. Additionally, high σ should have facilitated the emergence of broader virus host-range, which in nature has been shown to associate with longer interaction periods with hosts [30], potentially making these viruses strong limiting factors for their hosts. Emerging diversity within host and virus specific traits, such as host-range evolution in our model is a subject for a follow-up study. Preliminary results of diversity based on Shannon entropy in distributions of nutrient affinity, virus adsorption coefficient, and memory gene shown in the Supporting Information section indicate increasing diversity over the course of the simulations.

Overall, population dynamics are remarkably different for trade-off based vs random interaction functions. Long-term coexistence of both hosts and viruses is facilitated when the compatibility function is trade-off based (Fig 6a–6h), in which case the viruses persist the longest in closed systems (Fig 6e and 6g). Spatial sturcture has been discussed as mechanism facilitating long-term co-existence of bacteria and viruses [31], for example by rendering low virulent viruses higher fitness [32]. These findings might explain why eventually, even with compatibility trade-offs, virus populations die out in our well-mixed setting.

The dominance of variation in host traits over virus traits in regulating virus-host population dynamics described in the previous paragraph aligns with the observation that trade-offs involving host traits (ie, compatibility trade-off, which links host traits to virus traits, rather than virulence trade-off that is purely based on viral traits) appeared to have the most pronounced effects effects on the virus-host population dynamics. Indeed, further support for this is provided by our long-term simulations with and without compatibility trade-offs (Fig 6, top half vs bottom half), where trade-off based vs random compatiblity functions result in most distinct long-term dynamics, contrasting trade-off based vs random virulence runs in which difference are only pronounced in open systems (Fig 6, left half vs right half). Interestingly, long-term co-existence also most easily emerge when compatibility trade-offs are expressed. This is true even when viral pressure is high early on in the infection cycle in closed systems. We thus postulate that compatibility trade-offs facilitate virus-host coexistence observed in nature. Interestingly, observations of narrow host range viruses interacting with their hosts at lower total host abundance compared to broad-host range viruses [30] suggest that low host-range interactions are more efficient than those of high host-range viruses, which could be an expression of such a compatibility trade-off playing out in nature. Ranking traits and trade-offs involved in virus-host interaction according to their impacts on evolutionary outcomes is, however, challenging and remains to be studied further.

Our analysis on co-existence of host and viruses supports the hypothesis that trade-offs associated with competitive abilities of hosts and traits influencing the efficiency of viral infection are critical to shaping realistic long-term population dynamics and co-existence in virus-host systems. Competition and defense trade-off being a fundamental phenomenon in nature [33-35], we conjecture that findings from this study might be useful in understanding long-term perseverence of parasites in epidemiological context. Besides, in terms of conditions for high virus population numbers, results might also be useful in pharmaceutical applications, in particular large scale phage production for phage therapy. As a next step, we intend to shed light on host- and virus-population internal trait diversification as a consequence of this trade-off, touching upon the phenomenon of microdiversity in natural microbial communities.

Experimental data.

Data on bacteria—phage co-infection experiments rendered in Fig 3 of the main text. The data file inlcudes timeseries of optical density measured at 600 nm, showing number of replicates, mean and standard deviation with 15 minutes interval for different initial virus to host ratios used.

(CSV)

Click here for additional data file.

Model structures and algorithmic procedures.

The pseudocode describes details of initialization, interactions, budgets and evolutionary dynamics implemented in the model, step by step.

(PDF)

Click here for additional data file.

Model’s sensitivity to physiological parameters.

Host (a-h) and virus dynamics for various resource affinities for hosts (α), adsorption coefficients for viruses (β), maximum growth rates for hosts (μ) and burst size for viruses (κ). Plotted curves: ensemble averages from 100 runs for T = 36 h, with standard deviations shown as shading.

(TIFF)

Click here for additional data file.

Evolutionary diversification of host and virus genotypes.

Diversity time-series based on Shannon entropy of ensemble averages from 100 runs in a trade-off based compatibiltiy and virulence scenario shown for a) nutrient affinity of host (g), b) adsorption coefficient of virus (g), and c) memory gene of virus (g). Small subpanels show abundance distribution for the genotypes at three distinct time points over the course of T = 720 h.

(TIFF)

Click here for additional data file.

21 Jan 2022

Dear Dr. Våge,

Thank you very much for submitting your manuscript "Individual-based model highlights the importance of trade-offs for virus-host population dynamics and long-term co-existence" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Amber M Smith

Associate Editor

PLOS Computational Biology

Ville Mustonen

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Pourhasanzade et al. present an interesting manuscript studying how bacteria-phage interactions mediate long-term coexistence of the bacteria and phage. The authors develop a discrete-time individual-based model that includes potential tradeoffs between bacterial growth rates and their ability to defend against infection by phage. They compare modeling results to those of simple infection experiments of E. Coli and a phage at different VHRs on plates.

The population dynamics of E Coli are shown in figure 3. For higher VHRs, the bacteria dynamics are interesting in that the bacteria initially increase then decrease (presumably due to the bacteria-phage interaction) before increasing yet again later on around hour 20. Under certain parameter regimes, the model successfully recapitulates the qualitative nature of these curves.

The individual-based simulations include bacterial growth, death, division, mutation, and virus-host interactions. Generally, I felt that the description of the model could be clarified; perhaps a simple schematic that illustrates the various steps of the model would help a reader understand each element of the model.

I was confused by a couple aspects of the model. I did not see the time step given (is it one hour?). A couple of issues: the parameter pi_h is the probability of mutation of host genotype—surely this must depend on the time step chosen. The parameters delta_h and delta_v are stated as mortality (or decay) rates with units per hour in table 1, but then in the text description (line 76 and line 89) these same parameters are treated as probabilities.

While the model is of course highly idealized, it is still quite complex and contains a large number of parameters (table 1). Many of these parameters are fixed at very specific values given by the study in reference 21. While some sensitivity analyses were performed for a few parameters, it is difficult to know whether model results are specific to this system, or whether results are more general. A broader exploration of parameter space may clarify this.

A major component of the model is the inclusion of bacterial evolution. When a cell divides, its genotype changes with some probability. These genotypes determine the compatibility of the virus and the bacteria through the function C_t. Yet it was difficult to see the connection to data, or to disentangle the effects of evolution from the innate structure of the model (the ecology) in determining population dynamics. Figure 3 shows E Coli population dynamics, but I did not see any sequencing or anything to show whether these dynamics are driven by evolution on the time scale of the experiment. It is difficult to know to what extent evolution is driving these dynamics in the experiments. Moreover, all results in figs 4-6 show population dynamics over time. It is difficult to know to what extent evolution is driving these dynamics in the model. Also, phage population dynamics were not measured, so these cannot be directly compared to the results of the model.

Minor issues:

I would say that Nh and Nv are state variables and not parameters as stated in Table 1.

Typo on line 289 “nutrientt” instead of “nutrient”

I thought the figures could be improved by adding labels to the legends so that the reader would know what is being varied without having to read the caption.

Many of the figures are means of many replicates. It might be nice to shade to show the amount of variation in replicates (standard deviation or similar).

Reviewer #2: See attachment.

Reviewer #3: This paper investigates the hypothesis proposed by others (who are cited in the introduction) that a trade-off between bacterial host competitiveness and host defense against viruses will result in a persistent population of both virus and host. A model is described that takes into account mutation rate, nutrient availability and virus infection efficiency and the results compared to those of an in vivo experiment performed in E. coli with a T-4 like phage. The inclusion of the in vivo experiment is a major strengths of this work in my opinion.

The model parameters are clearly described in tables. The code is not included but is stated to be available upon request. Figures consist of graphs of parameters over time either in the live experiment or the model runs.

Figures are generally clear, but there are ways they or their presentation could be improved.

• It is stated that the model is validated with laboratory data from a short-term infection experiment. However, the in vivo experimental result (Figure 3) and model result (Figure 4a) for host abundance differ considerably at the lowest VHRs and this does not seem to be addressed in the text.

• Methods state that multiple repetitions were done, but no indication of variation between runs is indicated on the graphs and no raw data are provided.

The authors state clearly that their model is “highly idealized, focusing on specific trait-based mechanisms.” However I found myself disagreeing with two assumptions that the investigation is based on, and this lessened my enthusiasm for the work.

Assumption 1: Evolution (genetic diversification) of bacteria requires the long-term co-existence of bacteria and phages.

While it is true that predator-prey relationships are a source of selective pressure that leads to genetic change over time, it is not the only source of selective pressure and therefore not necessary for evolution. Competition for nutrients, as included in the model, are just one of many other selective pressures. The fact that viruses exist and affect population dynamics of their hosts, to me, is not evidence that they are necessary for a healthy ecosystem.

Assumption 2: There is always a trade-off between the (bacterial) host’s growth (competitiveness) and its abilities to defend against the virus. i.e. an increase in virus-resistance host will always result in slower reproductive rate.

If taken to the mathematical extreme, doesn’t this assumption imply that a host with 100% virus resistance should be completely non-competitive with regard to growth? But absolute resistance to phage is known to exist in viable bacterial strains, which means in a single virus-host pair scenario, the virus would go extinct. Doesn’t this disprove the single host-virus pair model?

minor typo: virues is missing an "r" in the abstract

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No: I did not see that the code for the simulations and the raw data from the laboratory experiment are made available.

Reviewer #2: No: Code is available upon request, but could probably be made available more immediately via SI or online repository.

Reviewer #3: No: The data points behind the plotted means and variance measures were not provided and there was no indication of a publicly accessible repository of the information or code.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: PLOSCB.pdf

Click here for additional data file.

23 Mar 2022

Submitted filename: CoverLetter_ResponseToReviewers.pdf

Click here for additional data file.

25 Apr 2022

Dear Dr. Våge,

Thank you very much for submitting your manuscript "Individual-based model highlights the importance of trade-offs for virus-host population dynamics and long-term co-existence" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Amber M Smith

Associate Editor

PLOS Computational Biology

Ville Mustonen

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: I thank the authors for their detailed responses. Most issues were adequately addressed. The authors now include supplementary pseudocode, a GitHub link to their code, and a more detailed model description, all of which make it far easier for a reader to follow their methodology. The figures are also greatly improved. One minor remains issue with the figures: Fig 3-5 are missing an explanation as to the meaning of the shading.

A key limitation that still exists is that phage population dynamics are not measured so that they cannot be directly compared to model results. In lines 284-287 of the discussion, this is touted as an advantage of the model. But I would argue that it is only an advantage if the viral population dynamics are independently verified with a rigorous comparison with data. Before such verification, the phage population dynamics seem speculative. In my view, such measurements are not needed for publication of this manuscript as the authors already state the challenges associated with obtaining such measurements. However, I do think that this limitation should be acknowledged.

Additional minor issue: Figure S2 does not appear to be referred to in the main text.

Reviewer #2: The authors seem to have addressed the majority of reviewer questions and concerns. While we might want a few more answers now, this paper seems to lay the groundwork for followup work.

I am satisfied with the revision and especially updates to the model description, tables/notation, and figures.

** I missed this the first time, but on line 145 a dot/bullet symbol is used in some chemistry notation. Not sure if this is standard.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

11 May 2022

Submitted filename: Response to Reviewers - Round 2.pdf

Click here for additional data file.

17 May 2022

Dear Dr. Våge,

We are pleased to inform you that your manuscript 'Individual-based model highlights the importance of trade-offs for virus-host population dynamics and long-term co-existence' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Amber M Smith

Associate Editor

PLOS Computational Biology

Ville Mustonen

Deputy Editor

PLOS Computational Biology

***********************************************************

2 Jun 2022

PCOMPBIOL-D-21-02137R2

Individual-based model highlights the importance of trade-offs for virus-host population dynamics and long-term co-existence

Dear Dr Våge,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol