Modelling the spatiotemporal complexity of interactions between pathogenic bacteria and a phage with a temperature-dependent life cycle switch

Model equations

We introduce a spatiotemporal model to describe the bacteria-phage interaction in the upper part of the soil with the depth H (we consider (H=1) m) in a typical agricultural field. Here we consider a 1D model where all abiotic and biotic components depend on time t and vertical coordinate h. The biotic component of the model consists of 4 compartments: phage-free bacteria (S) susceptible to infection by the phage, bacteria infected by the phage in its lysogenic ((I_1)) and lytic ((I_2)) states, and free phages (P). The total density of the host bacterial populations N is defined as (N = S + I_1 + I_2). The schematic diagram illustrating bacteria-phage interactions is similar to that of Egilmez and co-authors¹⁶. The local species interactions are described based on the classical modelling approach^6,19. Our spatiotemporal model is of reaction-diffusion type and is described by the following equations

In the above model, we parameterise the growth of susceptible bacteria via a standard logistic growth function⁶, where (alpha) is the maximal per capita growth rate and C is the carrying capacity of the environment; we assume that C(h) varies with depth. Infection of S by phages P at low temperatures results in lysogeny which is described by a mass action term (K_s S(t,h) P(t,h)). The growth of lysogenic bacteria (I_1) is described by a logistic function as in the case of S; however, with a different maximal growth rate ({overline{alpha }} (T)) as detailed in the next subsection. At high temperatures, the transition from the lysogenic to the lytic cycle of infection occurs: this is described by the term (lambda _1 (T) I_1(t,h)). Infection by the phage via the lytic cycle is modelled by the term (K_2 (T)S(t)P(t)). The death rate of infected bacteria due to lysis is modelled by (lambda _2 (T) I_2). The lysis of a bacterium results in the release of b new phages, the the burst size⁶. In the equation for P, KN(t)P(t) stands for the loss of phage due to binding to any type of bacteria (for simplicity, we assume that there is no saturation in binding). Finally, (mu P(t,h)) is the natural mortality or deactivation of phages.

According to this framework, the vertical displacement of the phage and bacteria are modelled by a diffusion term (first term in each equation), where (D_b) and (D_P) are the diffusion coefficients of bacteria and phage, respectively. The variation of the temperature T across the soil is described by the heat equation

where (D_h) is the diffusion coefficient of heat transfer (see more detail in the next section). Models (1)–(2) should be supplied with appropriate boundary conditions. We assume that the model has the zero-flux boundary condition for all biotic components (bacteria and phage) at (h=0) and (h=H). For the temperature, we consider Dirichlet boundary conditions such that (T(t,0)= T_s (t)) and (T(t,H)= T_H), where (T_s (t)) is the surface temperature and (T_H) is a constant temperature in deeper soil layers.

Parameterisation of equation terms

Next we describe the functional forms of the dependence of model parameters on the temperature and the depth. Following the previous study¹⁶, we assume that the maximal bacterial growth rates (alpha (T)) and ({overline{alpha }}(T)) are described by

where (T_0=38.2 ^circ text {C}) is the optimal temperature; (T_1=34.8 ^circ text {C}) is the temperature corresponding to the switch between the lytic and the lysogenic cycles; (alpha _{text{max}}=23 text {day}^{-1}) is the maximal possible growth, (sigma =9.1 ^circ text {C}) describes the decay of growth with temperature T^16,20.

In the equation for ({overline{alpha }}(T)), we assume that at a high temperature normal cell division of (I_1) stops since there is a transition to a lytic state in bacteria. In the soil bacteria grow anaerobically or microaerophillically, and the growth rates of B. pseudomallei under such conditions are yet to be studied. For simplicity they are assumed to be the same as under aerobic conditions. Realistic values of the above parameters are listed in Table 1. Note that in the model both (alpha (T)) and ({overline{alpha }}(T)) are in fact effective growth rates of the bacterial populations, i.e. they incorporate the replication of cells and as well as their mortality.

The overall adsorption rate of the phage K is estimated as (2 times 10^{-7} text {ml}^{-1} text { day}^{-1}) from Egilmez et al.¹⁶. The adsorption constants (K_1 (T)), (K_2 (T)) and the transition rate from lysogenic to lytic cycle (lambda _1(T)) depend on temperature as follows¹⁶:

where (K_S) is the maximal phage adsorption constant ((K_S=epsilon K) where (epsilon =0.3) is the adsorption efficiency) and (lambda _{1_text {max}}=23 text {day}^{-1}) is the maximal transition rate which is assumed to be equal to the maximal growth rate of the bacteria¹⁶. The lysis rate of bacteria (lambda _2=20 text {day}^{-1}) (depending on 50 min latency time¹³) and the burst size (b = 100) in the model are assumed to be constant¹⁶. The temperature dependence of (alpha (T)), ({overline{alpha }}(T)), (K_1 (T)), (K_2 (T)) and (lambda _1(T)) are shown in Fig. 2. The mortality rate of phages (mu) is high near the surface due to ultraviolet radiation, but the role of ultraviolet radiation becomes negligible starting from a depth of a few centimetres because sunlight cannot penetrate the soil. For the above reason, we can assume (mu =3 text {day}^{-1}) to be constant.

The carrying capacity C of the bacteria varies with the depth of the soil, according to empirical observations^21,22,23. This can be explained by the fact that the humus, oxygen, nitrogen contents, or/and water content in the soil generally decrease with depth²⁴. We use a combined approach to parameterise C(h) based on the available empirical data. We assume that in the absence of phages, the bacteria achieve numbers close to the carrying capacity at a given depth. Firstly, we parameterise the dependence of the overall bacterial load on depth in paddy soils in Southern Asia using the existing data²². Then we re-scale the obtained curve based on the available observations of B. pseudomallei at a depth (h=30 text {cm})^25,26. We approximate C(h) using the following simple Gaussian-type curve

where (C_text {surf}) gives the maximal number of bacteria near the surface (h), B determines how fast the bacterial abundance decreases with depth, (C_0) is background bacterial density which takes into account the fact that bacteria can survive even at large depths (e.g. (h=100 text {cm})). Based on our estimates (see supplementary material SM1 for more detail), we will use the following parameter values as defaults: (C_text {surf} = 1 times 10^6) (text {cell/ml}), (B=7.5 times 10^{-4}) (1/{text{cm}}^2), (C_0=10^4) (text {cell/ml}) . One can easily see that C(h) has a maximum at the surface and monotonically decreases with depth. We assume that the carrying capacity of the environment is not influenced by seasonal variation.

The coefficient (D_h) in the equation for the temperature distribution can be estimated as follows. Generally, (D_h) is related to (rho _s), (C_{rho s}) and (k_s) which are the bulk density, specific heat and thermal conductivity in soil, respectively, i.e. (D_h=k_s/(rho _s C_{rho s)}). We use the estimates for (rho _s), (C_{rho s}) and (k_s) from²⁷ which gives (rho _s=110.52 text {kg}/text{m}^3), (C_{rho s} = 1130) (text {J/kg K}) and (k_s = 0.0967) (text {W/m K}) and, for the diffusion coefficient (D_h=7.7 times 10^{-8}) (text {m}^2 text{s}^{-1}). The variation of (T_s)—the surface temperature—is obtained from the historical weather report for the surface¹⁶. The bottom boundary temperature (T_H) at (h=H=1 text {m}) is considered to be (22 ^circ text{C}). The initial value of the temperature distribution (T_s (0)) is assumed to be linear, but this assumption does not affect long-term temperature dynamics.

The paddy fields in which we model the bacteria-phage interactions are flooded lands, where the soil is either mud or muddy water. Many factors can affect vertical dispersal of bacteria and phages in such soil. For instance, rain water can carry bacteria and phage up or down in the soil, which can be mathematically modelled by adding an advection term; however, for simplicity we ignore such effects in this paper. We also assume the phage and bacteria vertical diffusion coefficients to be constant; however, it is rather hard to provide accurate estimates for (D_p) and (D_b). In water, the diffusion coefficient of bacteria and phages can be estimated as (3.6times 10 ^{-10} text {m}^2 text{s}^{-1}= 0.3 text {cm}^2 text{day}^{-1}) and (2.8 times 10^{-12} text {m}^2 text{s}^{-1}= 0.002 text {cm}^2 text{day}^{-1}), respectively²⁸, but the diffusivity in soil should be smaller than this. As such, these values should be considered as upper limits for (D_P) and (D_b), with the actual coefficients being orders of magnitude smaller. We undertook simulations with different combinations of diffusion coefficients in this range, and found that the patterns of vertical distribution do not largely depend on the diffusion coefficients provided (D_P< 10^{-3} text {cm}^2 text{day}^{-1}) and (D_b < 10^{-2} text {cm}^2 text{day}^{-1}), due to the strong external forcing of the system by temperature (see “Results” section for details).

In our numerical simulations, we use both explicit and implicit numerical schemes. We take a 0.1 cm spatial step size to get a proper resolution. We separately compute the heat equation to define T(t) with a smaller time resolution and then apply the temperature obtained to model bacteria-phage interactions for a larger time resolution (for example, (Delta t cong 7 times 10^{-5}) day). We compute the average densities of the species (both in terms of spatial and temporal averaging) using a numerical right Riemann sum. The accuracy of our numerical simulation was verified by reducing both time and space steps and comparing the results obtained. We use daily and seasonal variation of temperatures (for the period of 2013–2016) in the provinces of Nakhon Phanom and Sa Kaeo in Thailand to parameterise the model (http://www.worldweatheronline.com). The unit of the densities of bacteria and phages are cells/ml. The summary of model parameters as well their values are provided in Table 1.

Source: Ecology - nature.com