Pharmacodynamics of interspecies interactions in polymicrobial infections

Introduction

Polymicrobial infections (PMIs) involve the presence of multiple pathogens living in close proximity at the site of infection¹. It has been hypothesized that interactions between bacterial species, i.e., interspecies interactions, may contribute to the risk of treatment failure in PMIs^2,3,4. Interspecies interactions may alter pathogen sensitivity to antimicrobial drugs^5,6,7,8. Moreover, these interactions can impact evolution and establishment of antimicrobial-resistant mutants^{9,10,11,12,13}. The polymicrobial nature of various bacterial infections is increasingly recognized, in particular for respiratory tract infections in individuals with cystic fibrosis¹⁴, infections in patients with neutropenia¹⁵, blood stream infections¹⁶, chronic wound infections¹⁷, otitis media infections^3,13 and urinary tract infections^6,18. A better understanding of the complex dynamics of polymicrobial communities and their response to antimicrobials could help to develop enhanced therapies for patients with PMIs.

The impact of interspecies interactions on the bacterial population arises as a cumulative effect of numerous local interactions among individual bacteria. This impact depends on the strength and frequency of individual interactions. Interaction strength is inherent to the interaction’s mechanism. In contrast, interaction frequency can be shaped by various factors, including the degree of mixing between the two species and the distance at which interactions can occur. In spatially structured environments, local unmixing of species can occur as the result of growth and may limit the potential for interactions between individuals from different species at a short range^19,20. Examples of microbial populations living in such conditions are biofilm communities, very dense populations and those living in a movement-restricting environment, such as thick mucus. Conversely, planktonic bacterial populations composed of freely floating bacteria undergo continuous interspecies mixing which may prevent local unmixing and allow for more frequent interactions. Similarly, long-range interactions may be less sensitive to spatial effects, such as local unmixing, because of their potential to reach beyond a locally unmixed region. This underlines the importance of the interaction distance and the lifestyle of bacteria for their cumulative impact on the bacterial population.

Various interspecies interactions have been identified that influence sensitivity to antimicrobial drugs. These interactions encompass intricate and often poorly understood mechanisms, involving the exchange of metabolites^6,21, signaling molecules¹³, extracellular^5,22 and intracellular²³ drug inactivation, exoproduct secretion²⁴, environmental shifts²⁵, extracellular structures²⁶, virulence factors²⁷ or bacterial warfare²⁸. These interaction effects can be considered as a type of collective resistance, collective tolerance, or growth enhancement and impairment strategies⁹. Although interspecies interactions may differ in their mechanisms, they often lead to similar responses in the affected bacteria, narrowing the range of relevant interactions to be considered when studying bacterial drug responses. For instance, a decrease in drug concentrations surrounding the affected species could result from drug inactivation or physical barriers that impede drug distribution. Similarly, the replication rate of a bacterial species may be altered due to changes in nutrient availability within its environment or through a metabolic shift initiated by signals from other bacterial species. In such cases, although the experiences of the affected species may differ, if the outcome—such as a change in growth behavior—is the same, the impact on drug response may also be similar. This reasoning supports a generalized approach to categorizing interactions solely based on their effect on bacterial response to antimicrobials.

The pharmacological properties of antimicrobial drugs represent another key factor which may drive the impact of interspecies interactions on the bacterial response. Specifically, the interplay between the drug target and the interaction type becomes of interest. For example, antimicrobials that target replication and whose effectiveness is proportional to the replication rate of bacteria, like rifampicin, are likely to be influenced differently compared to antimicrobials whose effectiveness is not dependent on the replication rate, such as ciprofloxacin. Thus, considering drug characteristics alongside interspecies interaction types is important to develop a more comprehensive understanding of how interactions influence the response to specific drugs. This knowledge is necessary for the development of more effective therapeutic strategies that account for the complex dynamics of PMIs.

Studying PMIs using experimental approaches, both in vitro and in vivo, presents numerous challenges due to the complex nature of PMIs and the vast array of potential copathogens. While established methods and protocols for studying bacterial response to antimicrobials primarily focus on monomicrobial populations, novel techniques are being developed to study polymicrobial populations. Theoretical approaches offer the potential to gain insights into PMIs by integrating existing knowledge of eco-evolutionary dynamics and pharmacodynamics in silico²⁹. These theoretical models can inform future experimental designs and research efforts aimed at understanding how interspecies interactions impact the response of pathogens to treatment. Furthermore, in silico approaches have the potential to translate knowledge of how interspecies interactions influence pharmacodynamics into tailored treatment strategies.

In this study, we employ an in silico approach to explore the influence of interspecies interactions on the response of bacteria to antimicrobial drugs. To achieve this, we developed a model of an interactive multi-species bacterial community, combining an agent-based modeling (ABM) approach with ordinary differential equation (ODE) models. Specifically, we study two types of interspecies interactions that either transiently alter the growth or drug susceptibility of bacteria. We consider antimicrobial drugs with various pharmacological characteristics such as bactericidal, bacteriostatic, concentration-dependency and replication-dependency. Our model takes into account the heterogeneity within the bacterial population arising from these interactions, as well as the spatial distribution of individual bacteria.

We simulated unidirectional interaction scenarios between two pathogens, one of which is insensitive to drug treatment. Such a scenario can occur in polymicrobial infections where one pathogen is gram-negative and another gram-positive or when one of the pathogens is susceptible and another drug-resistant. Polymicrobial infections are usually treated with narrow-spectrum antimicrobial drugs, most of which are effective against either gram-negative or gram-positive pathogens. One well studied example of such a pathogen combination is Pseudomonas aeruginosa and Staphylococcus aureus which are common in lung infections of patients with cystic fibrosis and in chronic wound infections^30,31. Another example is Bacteroides fragilis in combination with Enterococcus faecalis or Clostridium perfringens which commonly co-occur in intra-abdominal infections³².

We systematically evaluated how interspecies interactions influence the bacterial pharmacodynamic response to antimicrobial treatment. Additionally, we investigated how spatial effects of interspecies interactions affect the pharmacodynamic response of bacterial populations, i.e., the maximal interaction distance and bacterial lifestyles with varying degrees of cell mixing, represented by the bacterial movement speed. Overall, this comprehensive study aimed to improve our understanding of how interspecies interactions shape bacterial responses to antimicrobial drugs.

Results

We studied the response of bacteria to antimicrobial drugs in the presence and absence of different interactions towards a focal pathogen. We simulated bacterial populations consisting of two bacterial species, a dynamic focal species and a constant secondary species. Interactions between them are affecting the focal pathogen’s drug sensitivity (i.e., MIC or maximal replication rate, called the interaction target). Interactions can increase or decrease the value of the interaction target, respectively called positive and negative interactions. To determine how interspecies interactions may influence response of the focal pathogen population to antimicrobial treatment, the focal pathogen population was analyzed at the time of maximal net growth.

The observed change in response of the bacterial populations in the presence of interspecies interactions is the cumulative response of all individual bacteria. Due to random spatial distributions of the bacteria on the grid, some individuals in a population may be exposed to more interactions than others, causing heterogeneity of drug sensitivity within the populations. Comparing the distributions of the interaction target characteristics (i.e., MIC or maximal replication rate; Fig. 1) in populations exposed to different drug concentrations showed that in unstructured populations the average change in a population is of similar magnitude at above-MIC drug concentrations and varies slightly only for concentrations close to the MIC (Fig. 2). The observed similarity is most likely due to the large number of bacteria unaffected by interactions and the selected time point.

Pharmacodynamics of interspecies interactions in polymicrobial infections — **Fig. 1: Interactions cause heterogeneity in the target population.**

**Fig. 2: The distribution of interaction target variability in an unstructured population is similar for all above-MIC concentrations.**

Interaction and drug type determine how pharmacodynamics are affected by interspecies interactions

We investigated the effect of interspecies interaction on pharmacodynamics of antimicrobials with five different pharmacological characteristics according to their concentration-response relationship (Table 1). The effect of interactions on pharmacodynamics, the interaction effect, was quantified separately for each tested drug concentration (see “Evaluation metrics” in the “Methods” section). We found that the presence of interspecies interactions can affect the antimicrobial concentration-response relationship in two ways: either the effect of interactions is present at all concentrations with similar effect size mimicking a vertical shift of the control curve, or it is only present at concentrations at which the drug effect is also sensitive to changes in drug concentration mimicking a horizontal shift (Fig. 3).

In the presence of interspecies interactions modulating growth, the pharmacological properties of the antimicrobials determined the type of effect observed in the concentration-effect relationship. The response to bactericidal drugs changed with a similar effect size at all tested concentrations compared to the control where no interactions are present (Fig. 3b, d). Interspecies interactions which negatively affect growth cause the pathogen to be more sensitive to the drug and lead to a decrease of the average maximal net growth rate. In contrast, positive growth interactions support the pathogen, causing it to be less sensitive, and lead to an increase of the average maximal net growth rate. With increasing strength of growth interactions, the observed interaction effect is larger for both negative and positive interactions. For bacteriostatic drugs, the response depends also on other pharmacological drug properties. The response to a proportional bacteriostatic drug is proportional to the drug effect only appearing in concentration ranges where the drug response is also sensitive to changes in the drug concentration (Fig. 3f). In contrast, the impact of growth interactions on the response to an additive bacteriostatic drug is similar as on bactericidal drugs (compare Fig. 3h to b, d).

Further, it is interesting to note that the change in population net growth rate caused by the negative growth interactions does not increase much with increasing interaction strength (Figs. 3b, d, f, h and 4b, d). The change in the net growth rate of the population is influenced by the interaction strength, the fraction of the population affected and how many interaction events occurred per bacteria. The impacts of the medium and strong interactions are very similar which is also reflected in the average change in the interaction target for the strong and medium negative growth interaction (Fig. 3e). This is caused by the choice of the interaction strength parameters and the natural limit for changes in replication rate and MIC, which is zero for both values. The natural limit of the replication rate is reached with a number of interactions than the natural limit of the MIC (compare Fig. 1g and e).

In contrast to growth interactions, MIC interactions impact bacterial response to all the drugs only in concentration ranges in which the drug effect is also sensitive to changes in drug concentration (Fig. 3a, c, e, g). The pharmacodynamics of concentration-dependent drugs is impacted in a similar fashion (Fig. 4). The simulated concentration-dependent drugs are formulated using the same equations differing from concentration-independent drugs only in the parameterization which explains the observed similarities (Table 1). However, unlike in the presence of concentration-independent drugs, the effect of the MIC interaction is also present at high concentrations of concentration-dependent drugs. For concentration-dependent drugs the drug effect has not reached its maximum at the upper limit of the tested concentration range and is thus still sensitive to changes in drug concentration.

**Fig. 4: Growth interactions impact the response to concentration-dependent drugs also at high concentrations.**

This similarity between the population-averaged changes in MIC and the replication rate for all tested concentrations in unstructured populations (Fig. 2) allows for prediction of the type of effect based on the relation of the drug and the interaction target. Studying the mathematical equations of each drug (Eqs. (4a)–(9b)) and the relation between the parameter affected by the interaction and the drug concentration matches with the two types of observed effects, the horizontal and the vertical shift (Fig. 3). The equations of drugs affected vertically have an additive relationship between the affected parameter, the maximal replication rate k_rep,max, and the drug concentration c which is in the drug effect E (compare Fig. 3b, d, h to Eqs. (4a), (4b), (4d) respectively). The equations of the drugs affected horizontally have a proportional relationship between the affected parameter, either the MIC or the k_rep,max, and the drug concentration (compare Fig. 3a, c, e, f, g to Eqs. (9a), (9a), (9b), (4c), (9b) respectively).

Interestingly, the different types of qualitative effects, vertical or horizontal shifts, show differences in the variability of the population net growth between iterations represented as shaded areas around the mean (Figs. 3 and 4). Differences in the magnitude of the interaction effect are caused by differences in the spatial distribution of bacteria on the grid. The standard deviation is most pronounced for strong positive interactions (Fig. 3a–e, g) and negative interactions that cause a large change from the control (Figs. 3a, c and 4c) apart from growth interactions in the presence of a bacteriostatic drug (Fig. 3f, h). For those scenarios which show a vertical shift (Fig. 3b, d) in the presence of a bactericidal drug, the variability increases with increasing drug concentration.

Variability of the mean value between iterations is a feature of heterogeneous bacterial populations. The control populations without interactions are homogeneous and therefore do not show much variability between iterations. We hypothesize that larger heterogeneity of the interaction target within a population is associated with larger variability of the mean population value between iterations. This explains that stronger interactions show larger variability. It further explains the small amount of variability of the negative growth interactions (Fig. 3f, h). Due to the natural limit of the replication rate and the parameter choice, the target heterogeneity caused by negative growth interactions is small compared to negative MIC interactions (Fig. 1). Additionally, we would expect high drug concentrations to more strongly select in favor of heterogeneity when positive interactions are present and against heterogeneity when negative interactions are present. This could explain the increase in variability with increasing drug concentrations observed for some scenarios (Fig. 3b, d).

In conclusion, the effect of interspecies interactions on drug pharmacodynamics can be predicted qualitatively based on the drug effect formulation best describing a certain drug because the average change caused by interactions in a population is similar at above-MIC concentrations.

Structured bacterial populations benefit more from positive interactions and are more resilient against negative interactions in the presence of bactericidal drugs

We investigated the impact of bacterial lifestyle (structured and unstructured) on the magnitude of the interaction effect and found that in the presence of bactericidal drugs, bacteria benefit more from positive interactions and are more resilient against negative interactions in structured populations than in unstructured ones (Fig. 5). While the bacterial movement in unstructured populations is undirected, the speed of movement is an important factor determining the amount of mixing of the two species and therefore the number of possible interactions between the species. In silico bacterial populations with fast moving bacteria represent well-mixed planktonic populations, called unstructured populations. Populations without movement represent very dense populations, biofilm populations or populations living in an environment that restricts movement and are called structured populations.

Fig. 5: The effect of positive interactions on the response to bactericidal drugs is enhanced in structured populations compared to unstructured ones, while the effect of negative interactions is weakened.

Comparing concentration-response curves of unstructured populations to concentration-response curves of structured populations showed differences in how strongly the interactions affect the bacterial sensitivity (Fig. 5a–d comparing the upper to the lower row). Visual analysis of the example for a bactericidal drug suggests that negative interactions of either type decrease the maximal net growth at certain concentrations less in structured conditions compared to unstructured ones. This shows that the bacterial population is less sensitive to the drug when the bacteria mix less because the impact of negative interactions is weaker. Further, we observed that positive interactions are affected in the opposite direction. In structured conditions positive interactions increase the maximal net growth at certain concentrations more than in unstructured conditions pointing towards an enhanced effect of the presence of positive interactions.

The magnitude of the interaction effect depends on the number of individual interactions occurring which is determined by the distribution of the two species on the grid. A spatial distribution that reduces the number of possible interactions will hamper the interaction’s overall effect; a spatial distribution that maximizes the number of interactions will enhance the effect of the interactions. In the beginning of each simulation, individuals from each species are placed randomly on the grid. This distribution then changes due to selective pressure of the drug and, in the case of unstructured populations, the movement of bacteria.

Quantitative comparison of the interaction effect in structured and unstructured conditions showed that structured conditions and the presence of a bactericidal drug enables the species distribution to change to the benefit of the bacteria over time by enhancing the effect of positive and weakening the effect of negative interactions (Fig. 5e). To quantitatively compare the differences between concentration-response curves in structured and unstructured conditions, the magnitude of the interaction effect relative to the control, the relative interaction effect (ρ), was determined for each concentration. It describes the percentual change in response of bacteria to the drug in the presence of interactions versus in their absence. Relative interaction effects in structured and unstructured conditions were then compared to each other using a log fold change. Whereby a positive log fold change indicates that the interaction effect is stronger in structured conditions compared to unstructured ones, and a negative log fold change indicates that the effect is weaker.

The response to bactericidal drugs is affected by the bacterial lifestyle (structured or unstructured) for both positive and negative MIC and growth interactions (Fig. 5e). The response to bacteriostatic drugs is unaffected except for the additive bacteriostatic drug. There are only minor differences between the proportional and additive drugs in the different lifestyle conditions except for the negative growth interaction for the additive bacteriostatic drug. Interestingly, a structured lifestyle is more beneficial for survival of the bacterial population in the presence of any type of interaction and any bactericidal drug; it enhances the effect of positive interactions and weakens the effect of negative interactions confirming the initial observation from Fig. 5a–d. Further, the concentration at which the largest difference was observed shifts to a higher concentration with increasing strengths of positive interactions and to a lower concentration with increasing strengths of negative interactions. However, in contrast to expectation, the log fold change for stronger interactions is not bigger than for weak interactions.

The differences observed for bactericidal and bacteriostatic drug can be explained by the different types of selection pressures they exhibit. When bactericidal drugs are present, bacteria that are in an unfavorable position, either receiving too few positive or too many negative interactions, are killed at a higher rate than those in more favorable positions. This creates a more optimal distribution of the species on the spatial grid at the end of a timestep than at the beginning. In structured populations, this distribution is conserved and in the next timestep few bacteria are likely to die because many of them are already in favorable positions overall benefiting the population by either maximizing or minimizing the number of individual interactions in the timestep. Bacteriostatic drugs however do only inhibit growth but not kill bacteria hence they have a much weaker selective effect. Therefore, bacteria in unfavorable positions will still survive and suffer again in the next timestep, cannot grow and then eventually die from natural causes in the presence of bacteriostatic drugs.

Long-range interactions can be more harmful to bacteria because they cannot benefit as much from the spatial rearrangement

We investigated the impact of the interaction distance on the magnitude of the interaction effect and found that long-range interactions can be more harmful to bacteria because they reduce the chance for an altered spatial distribution to appear (Fig. 6). The interaction distance is a crucial factor determining which species distribution is beneficial therefore influencing adaptation of the bacterial population and the magnitude of the interaction effect.

**Fig. 6: Long-range interactions are more harmful to bacteria because they cannot benefit as much from the spatial rearrangement to avoid or benefit from interactions.**

We considered a longer interaction distance of up to 200 μm because this leads to overlapping interaction groups compared to a short distance and is expected to affect spatial distribution. Importantly, the maximal distance of an interaction determines how many individuals are within a bacterium’s radius of interaction and how many individual interactions occur. A larger interaction radius, caused by long-range interactions, will include more bacteria than a smaller one. To investigate the effect of the interaction distance on the magnitude of the interaction effect, we corrected for the higher number of expected interactions in a perfectly mixed system when simulations are performed by adapting the interaction strengths (see “Interspecies interactions” in the “Methods” section). The observed differences of the magnitude of the interaction effect are thus only caused by the influence of the different interaction distances on the species distribution on the grid and not by differences in overall interaction strengths.

Comparing the effect of long-range interactions to short-range interactions when the bacteria are not moving reveals that a longer interaction distance can have an enhancing (indicated in blue), a weakening (in red) or no effect (in gray) depending on the type of interaction (Fig. 6a). We observed that at certain concentrations, positive long-range MIC interactions cause less of a sensitivity change than positive short-range ones at the same strength. Additionally, we found that negative long-range interactions are slightly stronger at medium and strong interaction strengths than short-range ones. This can be explained because a longer interaction distance decreases the chance for an interaction group with a species composition that optimizes the number of interactions to occur. The reason that a longer interaction distance hinders optimal spatial rearrangement is that the interaction groups are larger hence more average and that they overlap. In conclusion, in the presence of long-range interactions, bacteria have less chance to adjust to the presence of the interactions and drugs.

When bacteria are moving, the effect of positive long-range interaction on the response to bactericidal drugs is enhanced at some concentrations (illustrated in blue in Fig. 6b). The most affected concentration increases with increasing interaction strengths for the growth interaction. Further, a long interaction distance weakens the effect of weak negative interactions on the response to bactericidal drugs slightly (indicated in red), however it does not affect strong negative interactions. We also observed that medium and strong MIC interactions affect the response to bacteriostatic drugs less strongly when the interactions have a longer interaction distance.

Overall, the observed pattern in Fig. 6b shows similarities to both Figs. 5e and 6a which can be explained: Because of the longer interaction distance, the fraction of previous timestep neighbors which are also neighbors in the following timestep is increased compared to a short distance. This is also the case when comparing structured to unstructured conditions. In other words, when the interaction distance is longer, the bacterial population appears more static.

Discussion

We have developed a modeling framework which describes the dynamical response of multispecies bacterial communities to antimicrobials. Our analysis demonstrates that both the type of interaction and the type of drug determine how the pharmacodynamics of the focal pathogen are impacted and that this can be predicted qualitatively (Figs. 3 and 2). Furthermore, we found that changes in the bacterial lifestyle (structured or unstructured) and the interaction distance can influence the magnitude of the impact of interactions highlighting the importance of spatial organization of the species. The magnitude of the interaction impact is altered to the benefit of the bacterial population when the bacterial lifestyle, interaction distance and drug characteristics allow for adaptation of the bacterial population by changing its spatial distribution over time. In structured conditions, the bacterial population benefits more from beneficial interactions in the presence of bactericidal drugs and is more resilient against harmful interactions (Fig. 5). Similarly, when interactions have a longer interaction distance, beneficial interactions affecting the bacteria’s susceptibility are less beneficial while harmful interactions of any type are more harmful in the presence of any type of drug than when they have a short interaction distance (Fig. 6).

This works demonstrates that the effect of interactions can be predicted based on the mechanism of action of the drug in relation to the target of interaction (Fig. 3). This prediction is possible because the population-averaged changes of the MIC and the replication rate at the time of maximal net growth are similar for all tested concentrations (Fig. 2). This provides information of which type of antimicrobial may be able to benefit from or avoid the effects of interspecies interaction (Fig. 3). When an interaction affecting growth is dominant, the mechanism of action of a drug should be considered. At high concentrations, a bacteriostatic drug that acts proportionally to the rate of replication could be considered when the effect of interactions are to be avoided, while a bactericidal drug could be considered when the interactions are competitive and the treatment might benefit from the interaction effects. When interactions affecting resistance are dominant, we suggest that the concentration dependency of a drug should be considered when designing treatments (Fig. 4).

Our findings show that the bacterial lifestyle and the interaction distance can influence the magnitude of the impact of interactions highlighting the importance of spatial organization of the species. Previously, interspecies interactions have been shown to be a main driver of spatial organization of bacterial species^33,34. Further, the importance of spatial organization of bacterial species for community function, both in the absence^19,20 and in the presence of antimicrobial treatment^31,35,36,37 has been demonstrated. In the absence of antimicrobial drugs, it has been found that the amount of mixing between two species and the interaction distance are crucial for community function and growth in a multispecies cross-feeding community¹⁹. Conversely, in the presence of antimicrobials, it has been shown that antimicrobial exposure of a multispecies biofilm causes a spatial restructure of the community which promotes interspecies interactions protecting the bacteria against imipenem³⁶. Similarly, it has been shown that disturbances of the spatial organization in a biofilm can disrupt beneficial communication between bacterial species which in consequence decreases their expression of virulence factors³⁵. In comparing our findings with biofilm studies, it is essential to note that our model captures certain aspects of biofilm populations, particularly the restricted movement of individual bacteria within structured communities. However, it does not account for other characteristics that may affect the antibiotic sensitivity of the biofilm community. Overall, these findings highlight the potential of interventions targeting the spatial organization for antimicrobial treatment of PMIs and are in agreement with our findings.

Our analysis highlights the potential of additional therapeutics that target the bacterial lifestyle, interaction distance and the spatial organization of a multispecies community. More mixing in the population could decrease the effects of beneficial interactions between species and enhance harmful ones in the presence of bactericidal drugs (Fig. 5). Similarly, in structured microbial populations, an increase in the interaction distance could alter the magnitude of the effect of resistance interactions to the detriment of the bacteria because it reduces the chance that an optimal spatial arrangement develops (Fig. 6).

Our modeling framework for PMIs has an emphasis on interactions between individual bacteria and associated population heterogeneity. Our model offers a distinctive advantage by combining rule-based modeling techniques with deterministic ODE modeling. We used ODEs to describe the growth and kill rates of bacteria, while ABM techniques were employed to represent bacterial movement and interactions. The use of a stochastic approach such as the ABM facilitated the simulation of a heterogeneous population and interactions between individual bacteria in a straightforward manner. By incorporating ODEs, we limited the variability between iterations of the model solely to the variability arising from interspecies interactions which allows for easier interpretation of the simulation output and offered a clearer understanding of the impact of these interactions.

The developed model framework can be used to simulate pairwise or higher-order interaction networks within bacterial populations comprising multiple species in the presence of hypothetical antimicrobial drugs. Additionally, our model allows for the integration of pharmacokinetic models enabling the representation of temporal variations in drug concentrations. Moreover, the framework is readily extendable with other types of interactions, more detailed interaction mechanisms, multiple species, spatial heterogeneity of drug concentrations, such as drug gradients, and the action of immune cells. This adaptability opens up avenues for investigating a wide range of research questions related to polymicrobial populations and their response to antimicrobial treatment. One research area where our model holds promise is the study of pharmacological perturbations of microbiomes^38,39. By utilizing our framework, researchers can explore the effects of various pharmacological treatments on the dynamics and composition of microbiomes.

While our model’s design effectively supports the systematic analysis of interspecies interactions, it is challenging to compare our findings directly to experimental studies. This limitation arises partly from our deliberate decision to maintain a general approach and partly from the design of our in silico experiments. The scarcity of in vitro experiments focusing on PMIs restricts the experimental setup and measuring techniques available. Certain design elements in our simulations, which enhance result interpretability, cannot be replicated in an experiment. For instance, co-culturing two species under varied drug concentrations while maintaining a constant number of one bacterium is difficult to achieve experimentally. Our study aims to provide hypotheses, which can serve as a valuable starting point for future experimental research.

Designing experiments to validate the results derived in this study presents significant challenges. Key considerations include the selection of suitable bacterial strains based on their drug sensitivity, motility and specifically characterized interaction mechanisms. Further, the choice of measuring techniques is critical; the use of fluorescently labeled strains in coculture experiments at a single cell or population-level experiments will be essential. Depending on the experimental focus and questions, different approaches may be employed. Population-level experiments, aiming to assess outcomes such as drug response curves, may utilize co-culture assays to measure changes in total fluorescent signal of strains. Conversely, experiments examining spatial dynamics may employ imaging techniques to quantify in time spatial distribution and clustering of cell types in the presence of multiple drug concentrations. We believe that the model framework developed in this manuscript can provide valuable guidance for designing targeted experiments investigating the pharmacodynamics of polymicrobial infections, potentially reducing experimental effort.

The interpretation of our results is constrained by the simplified representation of interspecies interactions in our model and the choice of interaction parameters. Rather than considering the specific biochemical mechanisms underlying these interactions, our model focuses on their overall outcome. While most interaction parameters were set within a realistic range, determining the precise values of interaction strength and interaction distance from the literature was challenging. The mechanism of the interaction could provide important information for the parameterization of the interaction strength and the interaction distance. However, comprehensive information regarding the mechanism and quantitative properties of many interspecies interactions influencing drug sensitivity is currently lacking. As a result, we assigned them based on realistic population-level outcomes.

It is important to acknowledge that our study primarily explores unidirectional, pairwise interactions and assumes that one species is insensitive to drug treatment, even though complex higher-order, multi-directional interaction networks are likely to be present in PMIs. Our aim was to systematically analyze the impact of interspecies interactions on drug pharmacodynamics; thus we focused on the effect of interactions on one bacterial pathogen. It is important to point out that this assumption limits our ability to fully capture the intricate dynamics within PMIs. However, it is worth noting that our developed modeling framework is versatile and can be used to simulate bidirectional and higher-order interaction networks. We believe that our chosen approach, which unveils individual interaction dynamics, can contribute to shedding light on the intricate dynamics of higher-order interaction networks.

In summary, our study presents a comprehensive analysis of the impact of interspecies interactions on bacterial pharmacodynamics. The developed model is a versatile framework which can be easily adapted for a wide range of research questions related to polymicrobial populations and their response to pharmacological perturbations. Our findings offer insights into selecting antimicrobial drugs that can benefit from or avoid the effects of interspecies interactions and highlight the potential of interventions targeting spatial organization. Moreover, we provide hypotheses which can be a valuable starting point for future experimental research and may pave a way for the development of targeted treatment strategies that consider the complex dynamics of these infections, ultimately improving patient outcomes.

Methods

Model

We developed an integrated hybrid modeling framework that describes a microbial community consisting of multiple species exposed to antimicrobial drugs by combining an ABM approach and ODEs. An ABM approach was used to describe movement of bacteria and interactions between them, whereas ODEs were used to describe the growth of the bacteria and the effect of antimicrobial drugs. In the model, each bacterial cell in the population is represented by a computational agent, which we will refer to as bacteria in the remainder of this manuscript. The bacteria are positioned on patches of a two-dimensional grid with periodic boundary conditions representing an infection volume of approximately 1 nanolitre which can hold up to 1000 bacteria (Table 2 “Grid and movement”). The grid acts as a simplified representation of space providing information about which individual bacteria are in proximity to interact. During the simulation, bacteria can move on the grid in an undirected manner which represents passive diffusion.

Simulations of the model were performed in discrete, consecutive time steps. At every time step, each bacterium is activated and can perform four actions in a fixed order: it can move to a neighboring patch on the grid, interact with other bacteria, replicate or die as a result of natural causes or drug treatment (Fig. 7). Upon activation of a bacterium, each of the actions occurs with a predefined probability. First, the bacterium moves with probability p_move = n_steps. Second, the bacterium interacts with each bacterium of another species that is currently in its interaction distance with probability p_interact = 1. Third, the bacterium replicates or dies with probability ∣p_net∣ whereby positive values of p_net indicate a chance to replicate and negative values a chance to die. The value of p_net is individual to every bacterium and can be altered by interactions. It is derived from the deterministic population net growth rate k_net and the length of the simulated timesteps Δt according to:

$${p}_{{rm{net}}}={k}_{{rm{net}}}Delta t.$$

(1)

As a result, the dynamics of a bacterial population without interspecies interactions simulated in the ABM framework closely approximates deterministic dynamics described by ODEs (Eqs. (3)–(9b)).

Interspecies interactions

The model incorporates two types of interspecies interactions which represent interactions that affect individual bacteria’s sensitivity to antimicrobials that can occur when bacteria are in close proximity. The description of the interactions focuses on the outcome of the interactions, a change in different interaction targets T ∈ {k_rep,max, MIC}, but does not specify their mechanism of action. The growth interaction targets the maximal replication rate k_rep,max, while the MIC interaction targets the Minimum Inhibitory Concentration (MIC) of an individual bacterium. Each type of interaction can negatively or positively alter the interaction target of the affected bacterium, respectively called positive and negative interactions. In this project, the MIC was used as a measure of resistance of individual bacteria. The MIC of a population is the average of all individual bacteria’s MIC values which approximately corresponds to the experimentally used MIC measure signifying zero visible net growth after 24 h.

Interspecies interactions are characterized by their target T, their strength and the maximal distance at which they can take place. Interactions can be weak, medium or strong and occur at short or long range. During simulation, interactions occur between an activated bacterium and all bacteria of another species within interaction distance. Each interaction alters the initial value of the interaction target T₀ by a fixed fraction f according to its strength. While positive interactions cause an increase, negative interactions cause a decrease. When a bacterium is affected by multiple interactions in one timestep, the fractional changes are summed up for all N interactions to derive the altered value of the interaction target T_Δ of the bacterium according to:

$${T}_{Delta }={T}_{0}left(1+mathop{sum }limits_{i=0}^{N}fright),,text{for},,Tin {{k}_{{rm{rep}},{rm{max}}},{rm{MIC}}}.$$

(2)

The effect of an interaction is only present in the timestep in which the interaction occurs and is not inherited by daughter cells upon replication.

The parameter values of fraction f of MIC interactions, which define the strength of an interaction, were chosen so that the resulting changes in MIC range up to a four-fold change for a strong interaction. Fraction values of the growth interactions were matched to the MIC interaction for each strength to ensure compatibility between these interaction types. The change in the maximal net growth rate of a population in the reference conditions is the same whether a growth or an MIC interaction is present. The reference conditions are a unstructured population (with fast movement) interacting at short-ranges in the presence of a proportional, bactericidal, concentration-independent drug at a certain concentration.

The interaction strengths of long-range interactions were corrected to ensure that the overall effect is equal to the effect of short-range interactions under reference conditions. Thus, the fraction value f of a long-range interaction and its maximal value were normalized to the higher number of interactions in a well-mixed population. As a result, the observed differences in the magnitude of the interaction effect when comparing interactions with different interaction distances under conditions other than the reference conditions are caused by the influence of the different interaction distances on the species distribution on the grid and not by differences in the number of interactions.

The interaction distances for short- and long-range interactions were defined based on the number of bacteria which form an interaction group (Table 2 “Interactions”). When a bacterial population is growing at capacity in our simulations, an average interaction group of a short-range interaction encompasses 6.8 bacteria and half at the start of each simulation. This is approximately in line with experimental measurements reporting interaction groups consisting of 3 to a few dozen bacteria for short-range interactions²⁰. It is worth mentioning that the interaction distance measured in the same set of experiments is not comparable because of differences in the assumed bacterial density compared to our model. In the experimental setup bacteria are fixated in a microfluidic device in a very dense distribution. However, our in silico model allows for movement of bacteria, and we have assumed a lower maximal bacterial density based on polymicrobial lung infections of patients with cystic fibrosis⁴⁰ (Table 2). Further, we assumed an interaction distance that is twice as long and encompasses 9 times as many bacteria could be considered long-range compared to short-range ones.

Bacterial growth and drug effect

The net growth dynamics of the bacterial population exposed to antimicrobials in the absence of interspecies interactions were described by ODEs that were integrated in the ABM framework and summarized below. The growth of a bacterial population B consisting of multiple species is described by:

$$frac{dB}{dt}={k}_{{rm{net}}}B.$$

(3)

with the net growth rate k_net

$${k}_{{rm{net}},x}=left{begin{array}{ll}{k}_{{rm{rep}},Delta }-{k}_{d0}(1+{E}_{c,p})quad &,{text{for}}, {text{x}},={,text{c,p},}quad (4{rm{a}})\ {k}_{{rm{rep}},Delta }-{k}_{d0}-{E}_{c,a}quad &,{text{for}}, {text{x}},={,text{c,a},}quad (4{rm{b}})\ {k}_{{rm{rep}},Delta }(1-{E}_{s,p})-{k}_{d0}quad &,{text{for}}, {text{x}},={,text{s,p},}quad (4{rm{c}})\ {k}_{{rm{rep}},Delta }-{E}_{s,a}-{k}_{d0}quad &,{text{for}}, {text{x}},={,text{s,a},}quad (4{rm{d}})end{array}right.$$

determined from the replication rate k_rep, the natural death rate k_d0 and the drug effects E_x ∈ {E_c,p, E_c,a, E_s,p, E_s,a} which are incorporated in four different ways representing various types of antimicrobials (Table 1). The effect of bactericidal antimicrobials which kill bacteria, is applied to the natural death rate either in a proportional or an additive fashion and which are denoted by E_c,p and E_c,a, respectively. Similarly, to describe the effect of bacteriostatic antimicrobials which inhibit bacterial growth, the effect is applied to the replication rate in a proportional or an additive fashion and respectively denoted by E_s,p and E_s,a. We chose these four semi-mechanistic equations (Eqs. (4a)–(4d)) because they have empirically been found to describe the drug pharmacodynamics of various antimicrobial drugs correctly^41,42.

The total number of bacteria is limited by the carrying capacity K reducing the replication rate:

$${k}_{{rm{rep}},0}={k}_{{rm{rep}},{rm{max}},0}left(1-frac{B}{K}right)$$

(5)

at high numbers of bacteria. The replication rate altered by growth interactions is denoted by k_rep, Δ and similarly described as:

$${k}_{{rm{rep}},Delta }={k}_{{rm{rep}},{rm{max}},Delta }left(1-frac{B}{K}right).$$

(6)

The relationship between the drug effect E_x and the drug concentration C is assumed to be sigmoidal:

$${E}_{x}(C)={E}_{{rm{max}},x}frac{{C}^{h}}{{C}^{h}+{EC5{0}_{x}}^{h}}$$

(7)

with the Hill factor h, the maximal effect E_max,x:

$${E}_{{rm{max}},x}=left{begin{array}{ll}frac{{k}_{{rm{rep}},0}-{G}_{{rm{min}}}}{{k}_{d0}}-1quad &,{text{for}}, {text{x}},={,text{c,p},}\ {k}_{{rm{rep}},0}-{k}_{d0}-{G}_{{rm{min}}}quad &,{text{for}}, {text{x}},={,text{c,a},}\ 1quad &,{text{for}}, {text{x}},={,text{s,p},}\ {k}_{{rm{rep}},0}quad &,{text{for}}, {text{x}},={,text{s,a},}end{array}right.$$

(8)

and the drug concentration for which the effect is half-maximal EC50_x:

$$EC5{0}_{x}^{h}=left{begin{array}{ll}{{rm{MIC}}_{Delta }}^{h}frac{{G}_{{rm{min}}}}{{k}_{d0}-{k}_{{rm{rep}},0}}quad &,{text{for}}, {text{x}},={,text{c,p/a},},quad (9{rm{a}})\ {{rm{MIC}}_{Delta }}^{h}frac{{k}_{d0}}{{k}_{{rm{rep}},0}-{k}_{d0}}quad &,{text{for}},{text{x}},={,text{s,p/a},}.quad (9{rm{b}})end{array}right.$$

The relations of E_max,x and EC50_x have been derived analytically, similar to ref. ⁴³, from assumptions about the net growth rate at very high:

$$mathop{lim }nolimits_{cto infty }(psi (c))=left{begin{array}{ll}{G}_{{rm{min}}}quad &,{text{for}}, {text{x}},={{rm{c}},{rm{p}}/{rm{a}}}\ -{k}_{d0}quad &,{text{for}}, {text{x}},={{rm{s}},{rm{p}}/{rm{a}}}end{array}right.$$

(10)

as well as at MIC concentration. It is assumed that the net growth is zero when concentration equals the MIC concentration:

$$c={rm{MIC}}:psi ({rm{MIC}})=0.$$

(11)

Whereby G_min denotes the minimal net growth rate reached at very high concentrations of bactericidal drugs (for x = {c,p/a}).

Simulation design

Bacterial populations consisting of two bacterial species which can interact unidirectionally were simulated in the presence of antimicrobial drugs (for simulation settings see Table 2). Antimicrobial mono-treatment was simulated at constant concentrations ranging from a concentration equal to the MIC to 6 times its value for each of the hypothetical drugs (Table 1). Fifty iterations of the simulation were performed for each type of interaction, interaction strength, interaction distance, movement speed, drug type and drug concentration.

Table 1 Simulation of antimicrobial drug treatment

Full size table

Simulations of the model with different parameters for the movement, representing bacterial lifestyle, and the interaction distance were performed. An unstructured lifestyle, which is indicated by movement of the bacteria at every timestep, represents a planktonic population. Planktonic bacteria are not attached to each other or any surface and are freely floating and moving due to diffusion. The model is also simulated without movement (n_steps = 0) which represents structured bacterial population with severely restricted movement such as biofilms or very dense populations. We chose two settings for the interaction distance which represent short- and long-range interactions.

Each simulation run was initialized by placing in total 500 bacteria with equal amount of each species randomly on the grid (Table 1). During a simulation, bacteria belonging to the species that initiate interactions are insensitive to drug treatment. As a result, the number of bacteria from this species stays constant, which allows us to systematically investigate how the focal bacterial species is affected by the interactions and the drug. At the end of every timestep t of the simulation run, the following data were recorded: the number of bacteria B of each species j (B_j), the probability to divide or die (p_net,i,t), the maximal replication rate (k_rep,max,Δ,i) and the MIC_i of every living bacterium i.

Evaluation metrics

The drug response of the focal pathogen population was measured by the maximal absolute value of the population net growth rate k_net during the runtime of one iteration. The maximal net growth rate of a bacterial population is a commonly used measure of the antimicrobial effect because it is a direct measure of fitness⁴⁴. It was determined from the average stochastic net growth rate p_net of every bacterium i at time t as follows:

$${k}_{{rm{net}},{rm{max}}}:= mathop{max }limits_{t}({k}_{{rm{net}}})=mathop{max }limits_{t}(overline{{p}_{{rm{net}},i,t}})$$

(12)

and is denoted by k_net,max. The drug response to different drug concentrations was illustrated for each simulated scenario in a concentration-response curve.

The effect of the interactions on the drug response of the focal pathogen was quantified for every drug concentration to compare the simulated scenarios to each other and called the relative interaction effect ρ:

$$rho =frac{{{k}_{{rm{net}},{rm{max}},Delta }}^{c = x}-{{k}_{{rm{net}},{rm{max}}}}^{c = x}}{{{k}_{{rm{net}},{rm{max}}}}^{c = 6}-{{k}_{{rm{net}},{rm{max}}}}^{c = 1}}.$$

(13)

Therefore, the difference of the maximal net growth rate in the presence k_net,max,Δ and absence k_net,max at concentration c = x was normalized to the maximal possible change in the absence of interactions.

Further, a log fold change of the relative interaction effect was determined as:

$${log }_{2}left(frac{{rho }_{1}}{{rho }_{2}}right)$$

(14)

and used to compare the magnitude of the interaction effect in different scenarios.

Table 2 Model parameters and simulation settings