Computational Model To Quantify the Growth of Antibiotic-Resistant Bacteria in Wastewater

ABSTRACT Although wastewater and sewage systems are known to be significant reservoirs of antibiotic-resistant bacterial populations and periodic outbreaks of drug-resistant infection, there is little quantitative understanding of the drivers behind resistant population growth in these settings. In order to fill this gap in quantitative understanding of the development of antibiotic-resistant infections in wastewater, we have developed a mathematical model synthesizing many known drivers of antibiotic resistance in these settings to help predict the growth of resistant populations in different environmental scenarios. A number of these drivers of drug-resistant infection outbreak, including antibiotic residue concentration, antibiotic interaction, chromosomal mutation, and horizontal gene transfer, have not previously been integrated into a single computational model. We validated the outputs of the model with quantitative studies conducted on the eVOLVER continuous culture platform. Our integrated model shows that low levels of antibiotic residues present in wastewater can lead to increased development of resistant populations and that the dominant mechanism of resistance acquisition in these populations is horizontal gene transfer rather than acquisition of chromosomal mutations. Additionally, we found that synergistic antibiotics at low concentrations lead to increased resistant population growth. These findings, consistent with recent experimental and field studies, provide new quantitative knowledge on the evolution of antibiotic-resistant bacterial reservoirs, and the model developed herein can be adapted for use as a prediction tool in public health policy making, particularly in low-income settings where water sanitation issues remain widespread and disease outbreaks continue to undermine public health efforts. IMPORTANCE The rate at which antimicrobial resistance (AMR) has developed and spread throughout the world has increased in recent years, and according to the Review on Antimicrobial Resistance in 2014, it is suggested that the current rate will lead to AMR-related deaths of several million people by 2050 (Review on Antimicrobial Resistance, Tackling a Crisis for the Health and Wealth of Nations, 2014). One major reservoir of resistant bacterial populations that has been linked to outbreaks of drug-resistant bacterial infections but is not well understood is in wastewater settings, where antibiotic pollution is often present. Using ordinary differential equations incorporating several known drivers of resistance in wastewater, we find that interactions between antibiotic residues and horizontal gene transfer significantly affect the growth of resistant bacterial reservoirs.

IMPORTANCE The rate at which antimicrobial resistance (AMR) has developed and spread throughout the world has increased in recent years, and according to the Review on Antimicrobial Resistance in 2014, it is suggested that the current rate will lead to AMRrelated deaths of several million people by 2050 (Review on Antimicrobial Resistance, Tackling a Crisis for the Health and Wealth of Nations, 2014). One major reservoir of resistant bacterial populations that has been linked to outbreaks of drug-resistant bacterial infections but is not well understood is in wastewater settings, where antibiotic pollution is often present. Using ordinary differential equations incorporating several known drivers of resistance in wastewater, we find that interactions between antibiotic residues and horizontal gene transfer significantly affect the growth of resistant bacterial reservoirs.
KEYWORDS antibiotic resistance, mathematical modeling, wastewater W astewater and sewage systems are a major reservoir of resistant bacterial populations due to the collection of antibiotic waste from humans and animals, inappropriate drug disposal, and effluent from drug manufacturers, hospitals, and agricultural/veterinary settings (1). In some cases, environmental concentrations can reach, or even exceed, MICs of certain antibiotics (2). This problem is particularly pertinent in low-and middle-income countries (LMICs) where cases of antibiotic-resistant infections have been rising and 70% of sewage produced is estimated to enter the environment untreated (1). For example, in 2016, an outbreak of extensively drug-resistant (XDR) typhoid cases emerged in in the southern Sindh province of Pakistan, and geospatial mapping revealed that the XDR Salmonella enterica serotype Typhi infections spread around sewage lines in the city of Hyderabad in Pakistan (3). Such outbreaks put the health and safety of surrounding populations at risk and, due to the communicable nature of these infections, increase risk for populations beyond the immediate vicinity of wastewater and sewage lines. Antibiotic pollution in wastewater and sewage systems results in a complex environment with many interacting antibiotic residues, and this pollution is a major driver of the antibiotic resistance that is leading to outbreaks in surrounding communities (4,5). One reason for this is that selective pressure from antibiotics in the environment is known to promote chromosomal resistance mutations (2,6,7). Furthermore, the interactions between different antibiotic residues may also affect resistance due to the synergistic and antagonistic effects changing the selective pressure on the bacteria (8). Additionally, the concentrations of antibiotic residues and disinfectants in sewage may be able to support and promote horizontal transfer of resistance genes among bacteria (9)(10)(11)(12). These mobile resistance genes have previously been measured at high levels in wastewater and sewage (2). Though AMR and its drivers have been studied in sewage and wastewater settings to some extent (2,12,13), one of the largest gaps in understanding the emergence of AMR within a sewage environment is the limited quantitative understanding of the effects and interactions of the many biological and environmental mechanisms at work (14). Quantitative understanding of resistant population growth in wastewater settings is critical in order to predict the development of large resistant populations in wastewater that may pose a risk to local populations of resistant infection outbreak. Furthermore, quantitative tools for understanding of resistance can be used to develop and model strategies for preventing the emergence of large resistant populations and therefore influence policy decisions. Therefore, there is a critical need for developing quantitative methods of probing AMR in wastewater environments. Our study is a step in filling this gap in knowledge.
Mathematical modeling has been an important tool in quantitatively studying AMR development at both an epidemiological and mechanistic level (6,8,12,(15)(16)(17). Epidemiological studies have included approaches that investigate bacterial population dynamics in a number of biological contexts, including biofilms (15). At the population level, much of AMR modeling deals with disease states in which resistance patterns are long established, with significant gaps in the study of resistance emergence in new pathogens and the role of environmental factors in the development of AMR (16). At the mechanistic level, modeling has been used as a tool to understand the individual roles of both mutational resistance and plasmid-mediated resistance (8,12). Models based on resistance developed due to chromosomal mutations have been used to study the role of synergistic and antagonistic antibiotic interactions on the emergence of resistant populations, finding that increased synergy increases the likelihood of resistance acquisition (8). Additionally, studies using modeling as a tool to understand horizontal gene transfer (HGT) as a driver of AMR have found HGT, specifically through the mechanism of conjugation, to be a significant mode of resistance acquisition for Escherichia coli in settings, including agricultural slurry (12). However, much of the work done recently has been focused on fitting data in the absence of antibiotic exposure (17). Understanding how bacteria evolve upon antibiotic exposure is important to develop strategies that prevent the emergence of resistance (6). Thus, there remains a need for quantitative modeling of both chromosomal mutation acquisition and HGT under selective pressure from antibiotics to fully understand AMR development in settings such as sewage and wastewater where antibiotic residues are often present. Here, we present a model of emergence and growth of drug-resistant bacteria in wastewater, integrating acquisition of resistance through both chromosomal mutations and HGT while also incorporating the effects of antibiotic interactions (Fig. 1).

RESULTS AND DISCUSSION
Model validation. In order to validate the developed model, experiments were conducted with E. coli in a simplified experimental system including only one antibiotic (rifampicin) and no plasmid-containing bacteria such that there was only chromosomal mutation as a mechanism for resistance. This experimental validation was done on the eVOLVER, a continuous culture platform which allows for automated, highly flexible and scalable microbial growth and lab evolution. eVOLVER allows for continuous flow conditions and independent, precise, and multiparameter control of growth conditions such as temperature and flow rate (18). This allowed us to precisely recreate the model parameters in vitro in a high-throughput and repeatable manner. For these eVOLVER studies, we set certain inflow and outflow rates, initial bacterial populations, and antibiotic concentrations as required inputs for the computational model and then measured susceptible and resistant bacterial populations over the course of the continuous evolution experiment. The outputs of the study can then be directly compared to the outputs of the model. The results of these studies are shown in Fig. 2a to c. The experiments showed that at 6 mg/liter and 8 mg/liter rifampicin (where the MIC of rifampicin is 25 mg/liter), the resistant population developed to a steady state around the same order of magnitude of the steady-state population of susceptible bacteria. This steady-state resistant population was larger for the 8 mg/liter rifampicin condition. At 12 mg/liter, the experiments showed that the resistant bacteria dominated the population, and the susceptible population dropped to zero. After the conclusion of the experiment, model equations and parameter values were adjusted to match the experimental behavior (see Fig. 10 and Table 3). Notably, the mutation rates were increased by several orders of magnitude to match the rapid in vitro development of resistance. Additionally, the experimental results led to the addition of an equation for a theorized semiresistant population with a lower level of resistance (MIC = 100 mg/liter) (S 1 ). This was added to match the experimental finding that resistance can coexist at a steady state with the susceptible population. We theorized that the observed susceptible population also contained a subpopulation of semi-resistant cells which were not detected by the 8Â MIC resistance cutoff rate of the experiment. These semiresistant cells would allow for a sustained susceptible population in conditions with constant antibiotic pressure.
After adjustments, the resulting model outputs matched the general behavior of the resistant and susceptible bacterial populations in response to the selected concentrations of rifampicin ( Fig. 2d to f). This established the ability of the model to match in vitro behavior allowing for the development of a model with predictive capabilities.
Following this initial study, the adjusted model was used to make predictions for a second experiment with different rifampicin concentrations ( Fig. 3a and b). This second experiment was again conducted using the eVOLVER with the same growth conditions, and the results of the experiment qualitatively matched the behavior predicted by the model (Fig. 3c and d). This predictive capability demonstrated the validity of the modeling assumptions and prompted further simulation studies using the model.
Model simulations. Initial model simulations were done to compare the outputs of our model with known scenarios in E. coli under low antibiotic pressure. In the scenario of no HGT, as expected, only mutational resistance was observed (Fig. 4a). Likewise, for a scenario in which there is no chromosomal mutation (chromosomal mutation rate = initial mutant population = 0), only resistance from HGT was observed and due to an assumed higher growth rate due to a lower fitness cost for bacteria with resistance-conferring plasmids than those with chromosomal resistance mutations (19), this resistance overtook the susceptible population at a higher rate (Fig. 4b).
Very low antibiotic residue concentration can promote resistant population growth. We then turned our attention to simulating realistic scenarios in wastewater and sewage. We studied the effects of changing antibiotic residue concentration on the development of antibiotic resistance (Fig. 5). Low subinhibitory concentrations (0.5 to 2 mg/ml), as would be present in wastewater settings (4,5), were studied. These concentrations are an order of magnitude lower than the MICs of each of the two simulated antibiotics, which have been assumed as 25 mg/ml based on our experimental results with rifampicin. At these low concentrations, increased antibiotic residue concentration increased the rate of resistance development and decreased the time to resistant population dominance (defined as time when R m 1R p . S). This is in agreement with several prior studies linking subtherapeutic antibiotic levels with both chromosomal resistance mutation development and HGT (6,7,(9)(10)(11)(12). Future studies will look at higher antibiotic concentrations and the point of inflection where this increase in antibiotic concentration begins to prevent resistance through bactericidal activity. Furthermore, HGT was observed to be the dominant mode of resistance, in agreement with patterns reported in studies of E. coli in other settings (20).
Increase in HGT rate increases resistance at low concentrations. Additionally, the effect of HGT was modeled by increasing the HGT rate, b (Fig. 6). Prior studies have indicated that HGT rate is a less significant driver of resistance frequencies (21). However, our model shows that at very low subinhibitory concentrations of antibiotic, increased HGT rate significantly decreases the time to resistant population domination. This result indicates that increasing HGT rate allows resistance to be acquired at very low selective pressures where there are infrequent chromosomal mutations. This result is important as it shows that decreasing antibiotic levels in wastewater to low levels is not sufficient in preventing the growth of resistant populations. At very low concentrations of antibiotics, the HGT rate of antibiotic resistance genes in wastewater can affect the proliferation of resistant populations in bacteria with high gene transfer rates. As the HGT rates of the multitude of bacterial species in wastewater is not a commonly monitored parameter, this may be a factor to be mindful of in wastewater surveillance, as these rates can have a significant effect on the evolution of resistant population reservoirs even at low antibiotic concentrations.
Effect of bacterial killing rate. Because wastewater settings can often have low concentrations of antibiotics from various polluting factors (4, 5), we probed the specific effects of antibiotic residues on the development of resistant populations. Individual antibiotics have different killing rates based on factors such as their modes of action. Thus, we investigated the effect of bacterial killing rate on resistant population growth (Fig. 7). Increasing the killing rate terms (d max,1 and d max,2 ) was observed to decrease the time it takes for the resistant population to dominate under subinhibitory concentrations. This result is in agreement with previous  studies showing that increased selective pressure from antibiotics at subinhibitory concentrations can increase resistance development (6,7).
Synergistic antibiotic interactions increase resistant population growth. The interaction between two antibiotics has previously been shown to affect resistance acquisition (22). Synergy is the interaction of multiple drugs to have a greater killing action than the sum of their parts. Antibiotic synergy has also been shown to increase the likelihood of resistance acquisition at subtherapeutic doses (8). However, the effects of antibiotic interaction on the growth of resistant populations in wastewater settings, where many antibiotic residues can be present, has not previously been observed or modeled. In order to fill this gap, we have probed the effects of synergistic effects between antibiotic residues on the system by varying the synergy parameter, syn. Results from our model show that increased synergy between different antibiotics  Computational Model of AMR in Wastewater resulted in a decrease in the time to resistant population elimination at suprainhibitory concentrations (Fig. 8a). This result may initially seem counterintuitive but is in fact consistent with previous studies. This is because the synergistic action between the antibiotics increases the killing action to the point where it is effective on the resistant populations at lower concentrations than for antagonistic pairs. Hence, the synergy between the two antibiotics allows for greater bactericidal activity at lower concentrations. However, at subinhibitory concentrations, increasing synergy between antibiotics decreased the time to resistant population dominance (Fig. 8b). This is in agreement with other models in which synergy between antibiotics was observed to increase resistance acquisition (8). Our model further shows that these effects of synergy are observable even with low levels of antibiotic, such as the antibiotic residue concentrations present in wastewater. While the effects of antibiotic interaction on resistance acquisition have been previously modeled, low subinhibitory antibiotic concentrations such as those found in wastewater have not been well studied and are critical for understanding resistance development in wastewater settings (17).
Conclusion. Though we have been able to draw a number of conclusions on the relative effect of a variety of factors affecting resistance development in wastewater, we note that our model does have limitations. First, our model is based on parameter values from literature, not all of which have been experimentally validated against wastewater conditions. Additionally, other phenomena related to microbial growth have not been included in this initial model. For example, the phenomenon of dormancy can allow persister cells within the susceptible population to survive in the presence of antibiotic and resume growth after antibiotic removal, providing a method for susceptible populations to survive higher doses of antibiotic and develop resistance at a later time point (23,24). Additionally, biofilm can also serve as reservoirs of resistance and have been observed to be enriched with antibiotic resistance genes downstream from wastewater effluent discharge points (25,26). Finally, our model takes into account only the stationary and death phases; however, the addition of the long-term survival phase, which has been observed in bacterial species, including E. coli, could allow for more precise understanding of population growth dynamics (27). The addition of these terms in future modeling studies may allow for better quantification of resistance development. However, despite these limitations, experimental validation demonstrated the ability of our model to qualitatively predict in vitro behavior of bacterial resistance development in response to multiple subinhibitory concentrations of rifampicin (Fig. 3). We have also assumed a simplified case where the two mechanisms of resistance acquisition that are modeled are independent, that as interaction between HGT and chromosomal mutation and the behavior of cells having both are unknown. This assumption can be tested through future experimental work. Our model also does not account for multiple bacterial species or a multitude of antibiotic residues, as would be present in a complex environment like wastewater. Interaction between the multiple bacterial species as well as quorum sensing within bacteria of the same species could have an effect on resistant population size and is a future area of expansion for this model. Additional improvements to the accuracy and robustness of the model could be made through parameter initialization from field data and experimental validation of model inputs and outputs under specific wastewater conditions using the highly tunable eVOLVER platform.
Despite these limitations, our model provides new, and important quantitative insight on the evolution of resistant bacterial reservoirs. It also provides an integrated framework to incorporate several aspects of resistance acquisition and growth previously lacking from models focused on AMR in wastewater settings. In terms of results, we have shown that the HGT rate can be a significant driver of resistant population growth at very low antibiotic concentrations (Fig. 5 and 6). This indicates that HGT rates of bacteria in wastewater may be important to monitor in addition to antibiotic residue concentration. We have also shown that synergy between the antibiotic residues present in wastewater can increase the rate of resistant population growth, even at the low concentrations present in wastewater (Fig. 8). Thus, antibiotic residues in wastewater may pose a greater risk than might be expected without taking these antibiotic interactions into consideration. This has important implications for determining acceptable antibiotic levels in wastewater after treatment, as determining levels without considering antibiotic interactions may lead to overestimating the permissible Computational Model of AMR in Wastewater level of antibiotic and allow for the proliferation of antibiotic-resistant bacterial populations. This model can be adapted for use as a prediction tool for public health policy makers and be used to predict resistant population emergence in different sewage and wastewater conditions. Additionally, it can be expanded to be used to model different resistant outbreak prevention strategies in sewage and wastewater treatment.

MATERIALS AND METHODS
Our mathematical model of the growth of antibiotic-resistant bacterial populations in wastewater builds on prior studies (8,12,28) and extends it to incorporate a variety of critical, but overlooked, input factors specific to the bacterial species, antibiotics, and environments of interest. The inputs can be broadly classified into bacterial parameters, environmental parameters, and antibiotic parameters. Bacterium-specific input factors include the growth rates of antibiotic-susceptible and -resistant strains, mutation rates, and rates of HGT. The antibiotic-specific inputs, such as bactericidal activity and degree of synergy, allow for the study of the effects of drug quality and antibiotic pollution on the development of resistance. Additionally, environmental inputs, including physical fluid inflow and outflow rates and antibiotic residue concentration, allow for the modeling of resistance development in a variety of settings of interest. These input parameters can be used to model an output of resistant bacterial populations over time, thus allowing for the prediction of resistant population development (Fig. 1).
The model consists of ordinary differential equations governing the concentration of two antibiotics (C 1 and C 2 ) over time as well as equations modeling the susceptible population (S) and populations resistant to both antibiotics 1 and 2 from chromosomal mutation (R m ) or from HGT, specifically a resistance-conferring plasmid (R p ) over time ( Fig. 9 and Table 1). Additionally, two populations of bacteria that are resistant through chromosomal mutation to each antibiotic individually are modeled (R 1 and R 2 ). Each antibiotic is modeled with terms for the antibiotic residue concentration in the environment (E) and the antibiotic clearance rate (k e ). Each bacterial population is modeled with terms for growth rate (a), which is limited by the carrying capacity (N max ). The bactericidal activities of the antibiotics are modeled by the term for killing rate in response to each antibiotic (d max,1 and d max,2 ). These killing terms are modified by the antibiotic concentration where the killing action is half its maximum value for either susceptibility ðC 50 S Þ or resistance to each antibiotic ðC 50 R;1 ; C 50 R;2 Þ. These killing rates and half-max concentration terms are assumed to be equal for each antibiotic. Additionally, the killing rates are modified by a synergy term (syn) with syn , 1 indicating antagonistic interaction and syn . 1 indicating synergistic interaction. Also included in the model are terms for bacterial inflow (g s , g Rm , g R1 , g R2 , and g Rp ) and outflow (k T ) rates (Table 2). These inflow and outflow rates are determined by the physical flow rates of the system of interest, and the outflow terms were validated experimentally as described below. As  susceptible and resistant bacterial inflow and outflow rates for wastewater settings are not known, the model simulations assume bacterial inflow and outflow rates of zero. However, these parameters are still included due to their relevance for the modeling of systems such as wastewater where bacterial inflow and outflow will affect resistant population development and they may be quantified in future field studies. Chromosomal mutation is modeled through terms for mutation rates under antibiotic pressure (m) which are concentration dependent. Parameters governing mutation rate and antibiotic interaction are based on a previously developed model of Michel et al. (8). HGT is modeled assuming the mechanism of plasmid conjugation and a HGT rate (b). HGT model terms and parameter values are adapted from a model of HGT in agricultural waste of Baker et al. (12). Other model parameters are based on experimentally derived parameters of E. coli in piglet studies (28,29). We note that the incorporation of all of these parameters, which have previously not been studied at once, into one single model can provide insights into their roles in the emergence and development of resistant populations. These equations were coded and solved in MATLAB (R2018a, Mathworks Inc.). Experimental validation of the model was done using the eVOLVER system (24). The experiment was initialized with inoculation of LB medium with E. coli MG1655 in static conditions at 37°C. Then, inflow and outflow of rifampicin-containing LB medium at three concentrations (6, 8, and 12 mg/liter) were started at a flow rate of 16 ml/h. Following the beginning of inflow/outflow, each culture condition was sampled daily, and the concentration of total bacteria was monitored through optical density (OD) measurements taken through the eVOLVER system and resistant bacteria were calculated through plating on selective LB agar containing 200 mg/liter rifampicin. This rifampicin concentration is eightfold greater than the experimentally determined MIC for rifampicin (25 mg/liter). After the conclusion of the experiment, model equations and parameter values were adjusted to match the experimental behavior ( Fig. 10 and Table 3).