Historical Contingency in Microbial Resilience to Hydrologic Perturbations

Development of reliable biogeochemical models requires a mechanistic consideration of microbial interactions with hydrology. Microbial response to and its recovery after hydrologic perturbations (i.e., resilience) is a critical component to understand in this regard, but generally difficult to predict because the impacts of future events can be dependent on the history of perturbations (i.e., historical contingency). Fundamental issues underlying this phenomenon include how microbial resilience to hydrologic perturbations is influenced by historical contingency and how their relationships vary depending on the characteristics of microbial functions. To answer these questions, we considered a simple microbial community composed of two species that redundantly consume a common substrate but specialize in producing distinct products and developed a continuous flow reactor model where the two species grow with trade-offs along the flow rate. Simulations of this model revealed that (1) the history of hydrologic perturbations can lead to the shifts in microbial populations, which consequently affect the community's functional dynamics, and (2) while historical contingency in resilience was consistently predicted for all microbial functions, it was more pronounced for specialized functions, compared to the redundant function. As a signature of historical contingency, our model also predicted the emergence of hysteresis in the transitions across conditions, a critical aspect that can affect transient formation of intermediate compounds in biogeochemistry. This work presents microbial growth traits and their functional redundancy or specialization as fundamental factors that control historical contingencies in resilience.

Development of reliable biogeochemical models requires a mechanistic consideration of microbial interactions with hydrology. Microbial response to and its recovery after hydrologic perturbations (i.e., resilience) is a critical component to understand in this regard, but generally difficult to predict because the impacts of future events can be dependent on the history of perturbations (i.e., historical contingency). Fundamental issues underlying this phenomenon include how microbial resilience to hydrologic perturbations is influenced by historical contingency and how their relationships vary depending on the characteristics of microbial functions. To answer these questions, we considered a simple microbial community composed of two species that redundantly consume a common substrate but specialize in producing distinct products and developed a continuous flow reactor model where the two species grow with trade-offs along the flow rate. Simulations of this model revealed that (1) the history of hydrologic perturbations can lead to the shifts in microbial populations, which consequently affect the community's functional dynamics, and (2) while historical contingency in resilience was consistently predicted for all microbial functions, it was more pronounced for specialized functions, compared to the redundant function. As a signature of historical contingency, our model also predicted the emergence of hysteresis in the transitions across conditions, a critical aspect that can affect transient formation of intermediate compounds in biogeochemistry. This work presents microbial growth traits and their functional redundancy or specialization as fundamental factors that control historical contingencies in resilience.

INTRODUCTION
Interactions between hydrologic and microbial processes control the cycling of dissolved organic matter and other chemical substrates (Covino et al., 2018) in biogeochemical hotspots such as river corridors (Singer et al., 2016;Graham et al., 2019). Predictive biogeochemical modeling therefore necessitates development of reliable models of the dynamics of microbial populations and their interactions with hydrology. For this purpose, it is important to understand how hydrologic inputs alter microbial diversity, abundances, interactions, spatial organization, and functions.
While specific microbial phenotypes of interest are assumed to be definable based on known gene functions and environmental conditions (Jansson and Hofmockel, 2018), their dynamics often exhibit non-intuitive, complex dynamics under the influence of the history of perturbations, a phenomenon known as historical contingency (also termed ecological memory or legacy effect) (Fukami, 2015;Ogle et al., 2015;Vass and Langenheder, 2017). Despite its critical impact, we do not know much about how microbial functional responses to hydrologic perturbations are controlled by historical contingencies and how to account for these effects in microbial models (Hawkes and Keitt, 2015).
Historical contingencies of environmental microbial systems have been frequently observed in diverse contexts (Bouskill et al., 2013;Peralta et al., 2013Peralta et al., , 2014. For example, microbial enzyme activity and respiration in soil that increased with moisture were shown to be constrained by previous climate such as precipitation history (Averill et al., 2016;Hawkes et al., 2017). As another example, soil systems showed an enhanced protection of plants against a pathogen through the induced presence of disease-suppressive microbes by pre-exposure to the same pathogen in the past, just like an adaptive human immune system (Raaijmakers and Mazzola, 2016). In a recent study, it was shown that different land use significantly affected the transfer of soil organic carbon and bacterial diversity and functions in surface water runoff (Le et al., 2020). Historical contingencies related to hydrologic perturbations were reported by Goldman et al. (2017) where the biogeochemical response (i.e., CO 2 flux) of parafluvial hyporheic zone sediments to re-inundation was contingent on hydrologic history. Goldman et al. (2017) interpreted lower and higher CO 2 fluxes in sediment samples as being caused by relative dominance of fungal and bacterial species, respectively. All examples here indicate the importance of accounting for the role of microbial species controlling historical contingencies for predictive biogeochemical modeling. While a current trend is to increasingly incorporate microbial physiology and processes into biogeochemical models (Wieder et al., 2013;Wang et al., 2017), accounting for historical contingencies still remains a challenge due to the lack of understanding of controlling aspects in microbial systems, severely limiting our ability to reliably predict biogeochemical function (Widder et al., 2016).
As a critical step toward addressing this challenge, we designed an in silico study to examine historical effects on microbial responses to hydrological perturbations. Inspired by the field observations of Goldman et al. (2017) that linked biogeochemical historical contingency with the shifts in microbial populations in a community, we formulated a mathematical model of a microbial consortium composed of two competing microbial species growing with different growth traits in a continuous flow reactor. In the model, one species grows more effectively in lower flow rates and the other species in higher flow conditions so that they hold microbial growth traits as observed in Goldman et al. (2017). We used this model to test the hypothesis that trade-offs in individual growth strategies underlie historical contingencies observed at the community level In this work, we explored how hydrologic history and population dynamics lead to historical contingency of resilience, whereby microbial responses to future perturbations are dependent on past perturbations (Holling, 1996;Allison and Martiny, 2008;Shade et al., 2012;Song et al., 2015;Martiny et al., 2017). Our model revealed that historical contingencies in resilience can be more pronounced for specialized functions associated with individual species, compared to community functions to which member species commonly contribute, indicating the distribution of growth traits among member species plays a controlling role. The model also predicted the sophisticated dynamics of microbial communities including hysteresis (as an outcome of historical contingency), and endogenous dynamics (i.e., the continued changes in microbial populations in the absence of environmental variation). Beyond highlighting the importance of representing the dynamics of individual species for predictive microbial community modeling, this work provides a specific model structure for parameterization in the context of resilience and historical contingency toward predictive ecological and biogeochemical modeling.

Design of the Study
Before providing the details on model development and specific metrics used for quantifying the effect of hydrologic perturbations on microbial function in the context of historical contingencies, we here provide a brief overview of the design of our study. Our major goal is to examine how the history of hydrologic perturbations will drive the composition and functions in microbial communities depending on the microbial growth traits and their trade-offs; how this effect will be shown differently for redundant vs. specialized functions among member species. For this purpose, we formulate a mathematical model of a continuous-flow stirred-tank (i.e., well-mixed) reactor where two species competitively grow on the shared nutrient (S) and produce distinct products (P 1 by species 1 and P 2 by species 2). This setting enables assessing the impact of past hydrologic perturbations on microbial functions (i.e., historical contingencies of resilience) based on the following two categories: (1) consumption of S, which represents a function commonly performed among two species and (2) production of P 1 and P 2 , which exemplify unshared functions uniquely associated with individual species. To account for growth traits observed in Goldman et al. (2017), we design the kinetics of microbial growth such that species 1 grows faster in lower flow rates, while species 2 grows faster in higher flow rates. This allows us to examine the effect of the trade-off in growth along the flow rate as a key mechanism for the resilience of community functions, particularly those redundantly shared among member species (i.e., the consumption of S). Using the resulting model, we assess resilience based the time-to-recovery after perturbations. For generalizable analysis across scales, we provide the final model equations in dimensionless form. In the following, we begin model description with theories for the co-existence of competing microbial species, which is a pre-requisite in building a microbial community model.

Mechanism for the Coexistence of Competing Species
Implementation of species co-existence in a community model requires a special consideration because simple kinetics will end up with the dominance by a fast grower due to the competitive exclusion principle, a classical ecological theory asserting that only one species can survive the competition (Hardin, 1960). As one of the simplest mechanisms enabling the coexistence of competing organisms in a continuous flow reactor setting considered in this work, we chose densitydependent growth (including flocculation or crowding effect) (De Leenheer et al., 2006;Haegeman and Rapaport, 2008), while several other plausible mechanisms are also available, including: spatial heterogeneity , externally imposed temporal fluctuations Smith, 1981), endogenously induced temporal fluctuations through a specific type of interactions (Hsu et al., 1978;Butler et al., 1983), dormancy (Lennon and Jones, 2011), nutrient storage (Revilla and Weissing, 2008), exchange of metabolic by-products (Hesseler et al., 2006), and time delay in the nutrient uptake and growth (Freedman et al., 1989). Densitydependent growth is a reasonable representation of the growth of microorganisms in natural environments that often form biofilms in sediments and soils where the growth of inside species is limited due to lower availability of nutrients along the depth, indicating that the overall specific growth can be suppressed by crowding as the density increases (Roughgarden, 1971). Even in the case where no biofilm is formed, the density-dependent growth mechanism is still valid in situations where species do not effectively distribute in space because cells surrounded by others will have the same difficulty in accessing nutrients as the local population density increases.

Microbial Growth Kinetics
Based on the density-dependent growth mechanism, we modeled specific growth rates of species i (µ i ), which are split into nutrient-and density-dependent terms, i.e., µ ′ i (s) and ρ i (x 1 , x 2 ) as follows: and As shown in Equation (3), we imposed environmental constraints (such as limited space and resources) on both organisms to limit the total population of the community below a given carrying capacity (C). Species 1 was chosen to have faster growth kinetics in this work, i.e., µ ′ 1 (s) > µ ′ 2 (s). Therefore, we additionally considered the exponential term in ρ 1 (x 1 , x 2 ) so that µ 1 is suppressed as the population density increases because otherwise co-existence becomes impossible due to the dominance of species 1.

Reaction Stoichiometries and Mass Balances
We considered two reactions (R 1 and R 2 ), respectively driven by species 1 and 2 (X 1 and X 2 ) such that they compete for the shared substrate (S) to produce distinct products (P 1 and P 2 ): where stoichiometric coefficients Y S i and Y P i (i = 1,2) denote how many quantities of S and P i are consumed or produced whenever a unit mass of biomass is produced through R i . Specific values of Y S i and Y P i depend on the chemical formulae of S, X i and P i , as well as their mass units (i.e., mole or gram). For simplicity, we set all stoichiometric coefficients with unity values i.e., Note that the resulting reaction stoichiometries above equate the production rates of P 1 and P 2 with the growth rates of X 1 and X 2 , respectively. To mimic microbial growth in river corridor sediments, we considered a continuous-flow well-mixed reactor where the substrate is fed through water flow into the sediment as illustrated in Figure 1A. A continuous flow reactor configuration has been frequently used in the literature as a model system for addressing ecological questions [e.g., Kraft et al. (2014)] and provides advantages over batch cultures in studying mixed populations (Veldkamp and Jannasch, 1972). Under well-mixed conditions, dynamic mass balances for microbial populations, substrate, and products can be written as follows: where x i denote the population size of X i , s and p i are the concentrations of S and P i , t is the time, µ i denotes the specific growth rate (1/time) of X i (mass/volume), k d is the specific death rate (1/time) fixed as constant, D is the dilution rate (1/time), and the subscript in represents the quantities associated with the inlet feed flow, respectively. As implied by the same notation of k d for species 1 and 2, we assumed no difference in the specific death rate between the two organisms.

Non-dimensionalization
For scale independence, we non-dimensionalized the equations as follows: where Definitions of dimensionless variables and parameters used in the above equations were provided in Table 1. The values of kinetic parameters k x , K x , k y , and K y in Equation (13) were derived from the flocculation model by Haegeman and Rapaport (2008). Nutrient concentration in the feed flow in the last term of Equation (10) was set to 1 due to the normalization with respect to itself.

Metrics for Resilience
We imposed hydrologic perturbation by fluctuating the dilution rate following a rectangle-shaped wave. Based on the concept of engineering resilience (Holling, 1996;Song et al., 2015), we used "time-to-recovery" as a measure of resilience, i.e., how quickly a chosen system variable or function recovers the pre-perturbed state, while other metrics could also be considered (Hillebrand et al., 2018). Key variables and functions of interest in this work include: microbial and chemical compositions (i.e., x, y, and z) and their associated functions such as consumption rate of z (i.e., r z ≡ −fx − gy) and production rates of p (i.e., r p ≡ fx) and q (i.e., r q ≡ gy). These rates represent redundant (r z ) and unique functions (r p and r q ) among member species.
In quantifying the resilience of a chosen variable or function among those mentioned above, we ensured the system to be exposed to a cyclic perturbation of dilution rate for a sufficient time period; when the system reached a sustained oscillation, we stopped imposing perturbation and started to measure time-torecovery. We used two different metrics for the time-to-recovery as defined below: • Hitting time is when the variable or function hits the preperturbed level for the first time. We measured this as the first time when the variable hits an upper or lower boundary of the band around the pre-perturbed level. The bandwidth was defined by the 2 and 5% values above and below the pre-perturbed level. • Settling time is the time required for the variable or function to finally reach and remain within a given band of the pre-perturbed level. This quantity needs to be distinguished from hitting time when the system recovery takes place in an oscillatory way so that it may touch the boundaries of the band more than one time. In this case, the first and last times that the nutrient profile crosses the boundaries of the band correspond to hitting and settling times, respectively. If there are no oscillations during the post-perturbation period, hitting and settling times are identical. In measuring the settling time, we considered 2 and 5% bandwidths as defined above.

RESULTS
Using a model of microbial community growing in a homogenous continuous flow reactor (Figure 1A), we test the following hypotheses: (H1) distinct growth traits and their trade-offs along the flow rate gradient ( Figure 1B) leads to the shifts in the community composition subject to hydrologic perturbations with different intervals, (H2) the resulting microbial communities dominated by species 1 or 2 will show distinct responses to future perturbations ( Figure 1C), and (H3) the extent of historical contingencies in microbial functions depend on their redundancy or specialization, i.e., whether a given functions is performed in common among multiple species or uniquely associated with specific organisms.

Trade-Off in Microbial Growth Traits in Steady State
We examined how state variables including species abundances (x and y) and substrate concentration (z) vary along the dilution rate (δ) (Figure 2A), which are dimensionless variables of x 1 , x 2 , s, and D (Methods). The values of x, y, and z in Figure 2A denote their steady-state values that change as a function of δ.
For very low dilution rates (i.e., δ ≤ 0.05), species 2 (y) could not survive the competition with species 1 (x). In this regime, x linearly increased with δ as the supply of substrate accordingly increased. For higher dilution rates above the threshold, two species coexisted with the opposite dependence of their growth rates on the dilution rate, i.e., x decreased with δ, while y increased, so that their population curves intersect at δ ≈ 0.1. This shows that species 1 is more competitive in lower flow rate conditions, while species 2 can grow faster in the higher flow rate regime, confirming that the parameter setting in our model leads to the expected microbial growth traits as indicated in Methods (also see Figure 1B) Similar trends were observed for community functions such as consumption rate of the shared substrate (r z ) and production rates of distinct products (r p and r q ) ( Figure 2B). As the conversion rate of z (r z ) includes the contribution from both species that competitively consume the shared substrate, the function r z showed milder changes over δ ≥ 0.05, compared to specialized functions (r p and r q ) that are uniquely associated with only one species. As will be discussed more in the follow-up section, the trade-off in growth traits of the two species provides a mechanistic basis for enhanced resilience of the community. We used the steady state profiles in Figure 2 as a guide map to identify distinct operating zones where the system shows qualitatively different dynamics.

Compositional and Functional Responses to Hydrologic Perturbations
We performed dynamic simulations to examine how microbial communities respond to hydrologic perturbations. Disturbances were imposed by periodically changing the dilution rate (δ) between pre-chosen nominal (i.e., pre-perturbed) (δ norm = 0.09) and inundated conditions (0.19) (following a rectangular shape in time). For illustration, we simulated compositional and functional responses to the two specific cycles of perturbations with the inundation intervals of 3.5 and 9.5 imposed over 0 ≤ τ ≤ 200 (Figure 3). Among the three state variables denoting microbial and chemical compositions (x, y, and z), z showed the most significant fluctuation in time because the substrate was being fed into the reactor while microorganisms were not. As species 1 is a faster grower at the nominal flow rate (δ = 0.09) and species 2 grows faster following inundation (δ = 0.19) (Figure 2), the abundances of species 1 and 2, respectively decreased and increased under periodic inundation (Figures 3A,B). That is, in the case of the shorter inundation interval (i.e., 3.5), species dominance was switched over from x to y (Figure 3A), while the original composition was recovered after periodic perturbation was stopped. This result supported our hypothesis on the effect of hydrological perturbations on microbial populations (H1). Similar trends were observed for functional responses including production rates of products (r p and r q ) and consumption rate of the shared substrate (r z ). Temporal fluctuations of r p and r q were more significant than those of their populations (x and y) because r p and r q are functions of substrate concentration as well as species abundances.
Even in this simple simulation setting, we found several interesting dynamics: (1) species 1 (x) showed inverse responses (i.e., the initial change in its abundance in a direction opposite to the final outcome) when periodic inundation was stopped at τ = 200, while such over-or under-shooting was suppressed for species 2 and substrate concentration; (2) the deviation of r z during periodic perturbation was the least significant among others (when evaluated based on the minimum values of fluctuation); and (3) the recovery of r z was the fastest among others. The result (1) indicates that species 1 with relatively fast dynamics and high momentum can show complex dynamics (such as inverse response). The results (2) and (3) partially supported the hypotheses H2 and H3 by showing the highest resilience for r z , which may be ascribable to the functional compensation from member species with growth trade-off. In the next section, these two hypotheses are more rigorously addressed based on the quantification of time-to-recovery in response to perturbations.

Historical Contingencies in Resilience
We performed a more thorough analysis of how resilience is influenced by historical contingencies. Based on the concept of engineering resilience, i.e., the rate of return to the preperturbed state (Holling, 1996;Allison and Martiny, 2008;Shade et al., 2012;Song et al., 2015;Martiny et al., 2017), we evaluated the system's resilience based on time-to-recovery such as hitting and settling times with 2 and 5% criteria as defined in Methods. For δ norm = 0.09 (as previously considered), we examined how state variables and their recovery times change as a function of inundation interval (Figures 4A-C). In support of the hypothesis H1, state variables including microbial (x, y) and chemical compositions (z) changed as a function of dilution rate. This change was significant in the regime where inundation intervals were small (i.e., when perturbation frequencies were high) (Figure 4A), while the time-to-compositional recovery did not show such a sharp increase over the corresponding regime ( Figure 4B). For x and z, the hitting time (solid line) was not always the same as the settling time (open circle), meaning that these variables first entered the area within 2 or 5% range of the nominal (i.e., pre-perturbed) values and later escaped them prior to ultimately returning to pre-perturbation conditions (Figures 4B,C). The time-to-recovery for z was shorter than those for x and y as indicated by the result in the previous section.
Many of these trends observed from compositional variables (x, y, and z) were reproduced for functional variables (r p , r q , and r z ) (Figures 4D-F). As a notable difference, however, r z showed suppressed changes compared to z due to the functional compensation from member species as explained earlier. Consequently, the time-to-functional recovery for r z not only was significantly lower than those for r x and r y , but also showed a milder dependence on inundation interval (Figures 4E,F). These results provide a clear support for H2 and H3.

Hysteresis in Transition Across Two Distinct Perturbation Cycles
As shown in Figures 5A,B, we simulated sustained oscillations of compositional (x, y, and z) and functional variables (r p , r q , FIGURE 4 | The impact of perturbation interval ( τ ) on species abundances, substrate concentration, microbial functions, and their time-to-recoveries: (A) microbial (x and y) and chemical compositions (z), (B) time-to-compositional recovery calculated by the 2% criterion, (C) time-to-compositional recovery calculated by the 5% criterion, (D) consumption rate of the shared substrate (r z ) and production rates of distinct products (r p and r q ), (E) time-to-functional recovery calculated by the 2% criterion, (F) time-to-functional recovery calculated by the 5% criterion. In (B,C,E, and F), solid lines and open circles denote hitting and settling times, respectively. and r z ) at two different inundation intervals (5 and 15) to see how the transitions occur between them. Consideration of these two intervals (wider than those considered the previous case) was to ensure the dramatic differences in the fluctuation patterns so that one can see transition trajectories more clearly. We initially imposed the periodic perturbation with the interval of 5 and when sustained oscillation was reached (i.e., at τ = 200), switched to the interval of 15 -this process was repeated afterwards. Phase diagrams in Figures 5C-E showed the existence of hysteresis in the transition between the two periodically perturbed conditions with distinct intervals. The transition between the two limit cycles formed by periodic perturbations with intervals of 5 and 15 occurred through different trajectories. The hysteresis was most effectively shown for the phase diagram of the two compositional variables (x and y) ( Figure 5C).

DISCUSSION
A growing body of data indicate critical effects of historical contingency in biogeochemistry and microbial ecology, but theoretical frameworks that can predict the sophisticated behavior of environmental systems are still lacking. In this study, we formulated a mathematical model of a minimal complexity to reveal how historical contingencies in resilience to hydrologic perturbations are driven by microbial growth traits and their trade-offs. Our model enabled us to reveal important aspects of resilience and its historical contingencies as summarized below.
Historical contingency in the resilience of microbial systems was significantly affected by the degree of functional redundancy or specialization. Redundant functions commonly performed by multiple species (such as the consumption rate of the shared substrate in this work) were more readily recovered from perturbations, while specialized functions associated with individual species (such as the production rates of distinct products by species 1 and 2) showed slower recovery. Consequently, historical contingencies in resilience were more pronounced for specialized functions partitioned only into specific species. The enhanced resilience for redundant functions performed by majority of species is ascribable to their compensating roles in contributing to a total biomass stock that is more robust to disturbances than individual species. Related to this, we also found that the compositional recovery of microbial populations was prolonged even after a certain community function (e.g., the substrate consumption rate) became stable. While this endogenous dynamics was discussed as an intrinsic characteristic of compositionally complex natural communities (Konopka et al., 2015;Song et al., 2015), our model demonstrated that even structurally simple communities can develop those subtle dynamics in the process of stabilizing microbial functions against perturbations. Our work therefore provides a minimal structure of model for studying a mechanistic linkage between endogenous dynamics, functional redundancy, resilience, and historical contingency.
As an outcome of historical contingency, we also observed hysteresis, whereby the system takes different routes in transitioning across conditions. Hysteresis of biological systems has been typically reported based on forward and backward trajectories of state variables between multiple steady states (Kim et al., 2012;Gedeon et al., 2018), but the hysteresis in this work is defined in a more dynamic context, i.e., as being characterized by distinct transient paths between two sustained oscillations maintained across different inundation intervals. As a signature of historical contingency, a predictive understanding of hysteresis is critically important because it can lead to the unexpected formation of intermediate chemical compounds during the transitions when microbial communities undergo the dynamic shifts in their composition and functions through different trajectories. The lack of understanding of hysteric behaviors in microbial communities therefore may cause a significant error in model predictions, particularly for those associated with specific organisms, the populations of which undergo a more significant dynamic variation than the total community biomass.
Our work offers a direction toward building more mechanistic biogeochemical models. Stronger historical contingencies of specialized functions as predicted by the model highlights the importance of representing distinct traits of individual species for predictive modeling -because the lumped description based on community-level properties such as the total biomass and/or carbon use efficiency cannot predict complex aspects of biogeochemical function such as resilience and historical contingency driven by the interplay between distinct microbial species. In this regard, model extension though linking conceptualized microbial functions with specific biogeochemical pathways would be an important next step. Two reaction pathways driven by species 1 and 2 considered in our model can be translated as specific biogeochemical pathways, for example, such as denitrification and anammox, i.e., once NO − 3 is converted to NO − 2 (as the first step of denitrification), the subsequent reaction paths bifurcated afterwards are driven by distinct microbial functional groups. Therefore, this type of model that has both redundant (i.e., conversion of NO − 3 to NO − 2 ) and specialized functions (i.e., conversion afterwards through two pathways) would be suitable for investigating functionspecific resilience and historical contingencies as demonstrated in our study. For this direction, trait-based modeling approaches that appropriately combine the concepts of functional guilds and functional enzymes should be considered as an appropriate modeling framework (Allison, 2012;Bouskill et al., 2012;Song et al., 2014Song et al., , 2017. We expect that accounting for individual growth rates, and species-associated specific biogeochemical traits will significantly improve predictive capabilities for multiscale hydrobiogeochemical models.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ hyunseobsong/history.

AUTHOR CONTRIBUTIONS
H-SS developed the modeling concept and performed simulations. H-SS drafted out the article, which was edited by JS and EG. All authors contributed to the design of the research, read and approved the final manuscript.