
Modelling the effects of climate change on the interaction between bacteria and phages with a temperature-dependent lifecycle switch
Infection of bacteria by the phage at high temperatures (\(T \ge 37 ^\circ C\))
Our experimental results of infection of B. thailandensis E264 by the phage AMP1 at high temperatures (\(T \ge 37 ^\circ C\)) are presented in Fig. 1, and also in Table S2 in Supplementary Material SM2. They reveal that AMP1 plaqued with essentially the same efficiency at the temperature range from \(37 \ ^\circ C\) to \(39 \ ^\circ C\); however, efficiency of plaquing sharply decreased at higher temperatures, dropping by more than 3 orders of magnitude at \(40 \ ^\circ C\). No phage plaques were observed at \(41 \ ^\circ C\), and the only sign of bacterial lysis, probably caused by the lysis from without mechanism, was noticeable when a neat phage stock sample was used. It is important to note that B. thailandensis was able to grow and form bacterial lawns on all plates. Thus, these results suggest that AMP1 cannot exhibit lytic behaviour at temperatures exceeding \(40 \ ^\circ C\). This gives us an estimate for the critical temperature \(T_{cr}\) (where infection stops) in Eq. (1) to be \(T_{cr}\approx 40 \ ^\circ C\).
From the experiment, it also follows that the switch from infection to a non-infection regime occurs suddenly within a narrow temperature interval of approximately \(1.5 \ ^\circ C\), which gives justification to model the cease of infection by phage at high temperatures by the sigmoid function \(F(T,T_{cr})=\frac{T_{cr}^m}{T_{cr}^m + T^m}\), with \(m=55\), simularly to the way of modelling the switch between the lysogenic and lytic infection cycles at the temperature \(T_0\approx 35 \ ^\circ C\) (see the study by Egilmez et al.11).

Plaquing efficiency of AMP1 phage on bacterial lawns incubated for 24 hours at various temperatures. Each row represents a different 10-fold dilution.
The temperature and UV index variation in Thailand

Observed and forecasted temperatures and UV index for the Nakhon Phanom province, Thailand. (A) Daily maximum and minimum temperatures from historical data for the period from 2009 to 2023; (B) Daily maximum and minimum temperatures for the period 2024 to 2044 forecasted using the best fitted SARIMA model. We assume a linear trend in the temperature raise to be \(1.5 \ ^\circ C\) over 20 years. (C) Daily average UV index from historical data for the period 2009 to 2023; (D) Daily average UV index for the period 2024 to 2044 forecasted using the best fitted SARIMA model. We assume the existence of a linear trend in the raise of UV index to be 3 units over 20 years.
Our analysis of the available meteorological data demonstrates a raise of both the mean temperature and mean UV index between 2009 to 2020 in the eight considered provinces in Thailand (see Fig. S1 from Supplementary Material SM3). This corresponds to the reported global land and ocean surface average temperature anomaly since 200037. This global anomaly ceased in years 2021 and 2022, with a significant drop in both the temperature and solar radiation.
Further, our investigation of the seasonal variation of temperature and the UV index shows that both quantities exhibit clear seasonal cycles, as shown in Fig. 2 (see panels A and C), constructed for the Nakhon Phanom province, as an illustrative example. One can also see that the temperature and the UV index also show pronounced variations on a daily basis. Figues also show that the peak of the temperature and UV index usually occurs in April and May. We should stress that the other considered provinces in Thailand also demonstrate similar patterns of variability in the temperature and UV index.
Finally, Fig. 2 (panels B and D) provides an example of a forecasted pattern of the temperature and UV index variation for the period 2024–2044. Note that in the shown forecast, we assume a linear trend in the raise of the UV index to be 3 units over 20 years. The provided forecast is for the Nakhon Phanom province; for the other considered provinces of Thailand we obtained similar results.
Modelling bacteria-phage dynamics during the period of 2009–2023
Using the mathematical model and the historic time series for environmental conditions (the temperature and the UV index), we performed numerical simulations of the bacterial and phage dynamics in Thailand for the period 2009 to 2023. Figure 3 shows the annual average values for the four model components \(S,P,I_1,I_2\) for the eight considered provinces of Thailand. Our simulations reveal pronounced inter-annual variability of the bacterial and phage densities over the considered period of 15 years. For most provinces, high values of phage-free bacteria S are found to occur for years 2013–2016 and also for 2018–2019; the lowest value of S is observed for year 2022. The mentioned inter-annual variability in S can be explained by the corresponding variation in the temperature, the total annual number of hot hours (with the temperature \(T > 35 \ ^\circ C\)) ‘HH’ and the UV index, which can be seen from the calculated Pearson’s correlation coefficients. The corresponding correlation matrix is provided in Supplementary Material SM4 (Fig. S2 and Table S3). Our estimates of pairwise correlation coefficients between the environmental components (T, UV, HH) and the bacterial density S are: \(\rho _{T,S}=0.71 \pm 0.08\); \(\rho _{UV,S}=0.90 \pm 0.03\) and \(\rho _{HH,S}=0.83 \pm 0.04\). In particular, the largest suppression of the bacterial density S in year 2022 can be explained by the drop in both the UV index and the temperature, signifying a more intensive control of bacteria by phages due to their reduced mortality, lower growth rates of the pathogen as well as the fact that at cooler temperatures (i.e. less amount of hot hours), infection by the phage of a bacterial cell results in lysogeny rather than in lysis. We should also stress that the inter-annual variation of \(S,P,I_1,I_2\) shown in Fig. 3 has a pronounced synchronisation across the considered Thailand provinces, which is due the synchronicity of the changes in the environmental factors over considered geographic region (cf. the curves in Fig. S1 from SM3).

Simulated dynamics of bacterial-phage interactions for years 2009–2023 for the considered eight provinces of Thailand (based on a non-spatial model). In each panel, the annually average values of species densities are plotted for different years. The names of provinces are indicated in the figure label. The unit of the densities of bacteria is cell/ml and phages is pfu/ml.
Typical patterns of the annual (seasonal) variation of bacterial-phage interactions are presented in Fig. 4, which is constructed for the temperature and UV data corresponding to the Nakhon Phanom province, as an example. In this figure, panels (A)–(C) provide dynamics for three selected years (years 2019, 2021, and 2022), whereas panel (D) shows the variation of species densities averaged over the total considered 15 year period. Figure 4 schematically shows a typical seasonal schedule of agricultural works in rice fields. In particular, the period of rice planting is denoted by the green rectangle, whereas the yellow rectangle represents the harvesting period. One can see that the model predicts a pronounced seasonal variation of the density S of susceptible bacteria as well as that of free phages P. High proliferation of S is mostly caused by a drop in the number of free phages P, which, in turn, is due to hot temperatures, where no infection of bacteria is possible (\(T>T_{cr}\)), and due to high phage mortality related to increased UV radiation. The annual variation of lysogenic bacteria \(I_1\) exhibits the inverse pattern as compared to that of susceptible bacteria S. This is explained by the competition between \(I_1\) and S for the resources, which is described in the by the common carrying capacity C. The density \(I_2\) of bacteria in a lytic stage is almost constant throughout the year, since the daily temperature variation allows for both lytic and lysogenic infection types each day, which results in a quasi-constant (on average) supply of \(I_2\). Major peaks of S, which present the most danger for humans, usually occur during the April-May period, which overlaps with the start of the agricultural works. Comparison of panels (A)–(C) demonstrates that the intensity of peaks of S largely varies from year to year, i.e. there exist particularly dangerous years (e.g. see panel (A)), where the density of the phage-free pathogen remain around the highest possible value, the carrying capacity C, for several weeks.
Our simulation also predicts high amplitude daily oscillations of species density. An illustrative example is shown in Fig. 5, constructed for the Nakhon Phanom province, and corresponding to the beginning of May, which is the start of the rice planting season. We present daily variation of \(S, I_1, I_2, P\) for May 1st, 3rd, 5th and 7th of year 2021. The model suggests that in the afternoon bacterial infections by phages are mostly lytic. On the contrary, in the morning, late evening and at night, phage infections mostly result in bacteria lysogenisation. Importantly, the susceptible phage-free bacteria S attains its peak in the evening (from 5 p.m. to 7 p.m.), and has its minimal density around noon. This makes the evening hours the most dangerous for field works. From the figure, one can also see a high variation of bacterial abundance on a day-to-day basis. This can be explained by the day-to-day variation of the ambient temperature and UV index. Similar patterns of daily variation of \(S, I_1, I_2, P\) are predicted by the model for the other provinces and years (the corresponding graphs are not shown here for brevity).

Modelled seasonal variation of the bacteria-phage system for the Nakhon Phanom province (a non-spatial model). Panels (A–C) provide the dynamics corresponding to years 2021, 2022 and 2019, respectively. These represent examples of typical years with a high, low and a moderate level of S, respectively. Panel (D) shows the dynamics averaged over the total period for 2009–2023. The green and yellow rectangles indicate typical seasons of agricultural work in rice fields: planting and harvesting, respectively. In each panel, the curves are plotted based on the moving average window technique with the window size of 30 days. The unit of the densities of bacteria is in cell/ml and phages is pfu/ml.

Simulated daily variation of bacterial and phage densities at the beginning of May 2021 (May 1st, 3rd, 5th, and 7th) in the Nakhon Phanom province (a non-spatial model). Panel (A) shows the dynamics of S and P, whereas panel (B) shows the dynamics of the infected bacteria \(I_1\) and \(I_2\). The unit of the densities of bacteria is cell/ml and phages is pfu/ml.
Forecasting bacteria-phage dynamics for the period 2024–2044
Figure 6 shows simulated model predictions of the bacterial and phage annual average numbers over the next 20 years, starting from year 2024. These predictions use our forecast of the temperature and UV index. We show the results obtained for the Nakhon Phanom province, as an example: our simulations for the other provinces provide similar patterns. The predictions are obtained for 5 scenarios of future trends in the UV index. In Fig. 6, we symbolically denote different trends in the UV index by UV(i), \(i=1,2,3,4,5\). This describes a range of alteration in the intensity of solar radiation from a light to a strong increase. For all considered scenarios, the temperature trend is kept to \(1.5 \ ^\circ C\) over 20 years.
One can conclude from Fig. 6 that for a light increase in UV radiation (scenario UV(1)), both the susceptible bacteria and phage densities show small growth. For a moderate increase in the UV index (scenario UV(2)), the density of phage is predicted to decrease, whereas that of bacteria shows pronounced growth. A stronger positive trend in solar radiation (scenarios UV(3), UV(4)) results in a substantial drop of P to the point where the phage goes extinct by the end of the considered time period. The density of susceptible bacteria eventually reaches its carrying capacity, i.e. the largest possible value, and the density of infected bacteria (both lysogenic and lytic) drops to zero. This is related to a pronounced increase in the mortality (deactivation) rate of phages by a very strong UV radiation to the point where an increase of susceptible bacterial host numbers would not compensate the loss due to phage mortality. Finally, for a very strong UV index (scenario UV(5)), the density of susceptible bacteria starts decreasing towards the end of the period of our forecast. This is related to the fact that too strong a level of UV becomes harmful for bacteria. Therefore, at the high levels of the UV index, pathogenic bacteria are regulated by solar radiation rather than by the phage.
Interestingly, the global raise in the temperature by \(1.5 \ ^\circ C\) over 20 years has much less impact on the bacteria-phage interactions, than the increase in the UV index. This can be seen, for example, from Fig. 6 for the scenario UV(1): the resultant increase in the density of S over the whole period due to the temperature raise for a fairly constant UV level is only about 20%. We also briefly investigated the dependence of our prediction on other model parameters. We found that an increase in the carrying capacity C, for example due to extensive use of fertilizers, would amplify the density of susceptible bacteria, resulting in increase of risk of disease acquisition. For high values of carrying capacity, an increase in the temperature and UV index results in the occurrence of outbreaks of high density of S throughout the year towards the end of the forecast period. The corresponding graphs can be found in Supplementary Material SM5. We also found that an increase in the burst size would delay the extinction of phages (we do not show the corresponding graphs for brevity).

Predicted outcomes of bacterial-phage interactions for years 2024–2044 for the Nakhon Phanom province (based on a non-spatial model). In each panel, the average annual values of species densities (\(S,P, I_1, I_2\)) are plotted for different years. The predictions are obtained for 5 different scenarios of the global trend in the UV index (denoted by UV(i), \(i=1,2,3,4,5\)), corresponding to an increase of the UV index over 20 years by 0.25, 3, 6, 6.5 and 8 dimensionless units, respectively. The unit of the densities of bacteria is cell/ml and for phages is pfu/ml.
Modelling effects of implementation of phage-killing agrochemicals.

Effects of the implementation of phage-killing agrochemicals on the bacteria-phage dynamics (using the non-spatial model). Free phages are removed from the system on September 1st. In each panel, we show the dynamics of the original (no phage killing) and the perturbed system (with phage killing) in the Nakhon Phanom province. Panels (A) and (C) are constructed for the year 2019, whereas panels (B) and (D) are obtained for the year 2021. The unit of the densities of bacteria is cell/ml and phages is pfu/ml.
Implementation of certain agrochemicals, in particular those with a high concentration of iron (II), can cause a severe mortality of phages21. In this section, we explore how implementation of phage-killing agrochemicals would affect the risk of infection by the pathogen under realistic fluctuation of weather conditions. We should stress that in this section, we model the effects of phage killing in the surface water of the rice field. We also briefly consider the effects of phage killing on bacteria-phage dynamics in soil.
In our model, we set the free phage density P to be zero on the day of release of chemicals. As the date of release, we select September 1st each year, and consider the period from 2009 to 2023. We compare the levels of S and P in the case of phage removal with the original system, where no phage-killing agrochemicals are applied. Examples of dynamics of S and P, in the case of removal of free phages, are shown in Fig. 7, constructed for the Nakhon Phanom province. In the figure, the left and right columns, correspond to years 2019 and 2021, respectively. The figure demonstrates two different possible outcomes of killing free phages. Namely, phage removal may have only a small impact on the system, as shown in panels (A) and (C): the species densities in the original and the perturbed systems are close to each other; the level of susceptible bacteria is generally low. Another scenario is demonstrated in panels (B) and (D) of Fig. 7, where removal of free phages results in a pronounced increase in the amount of phage-free bacteria. The level of S attains its maximal possible value, the carrying capacity. Importantly, for both above mentioned scenarios, the system demonstrates strong resilience: the initial perturbation of the system due to phage killing eventually becomes negligibly small, so the system returns to its original (i.e. non-perturbed dynamics) course of time. This also indicates the external forcing by the temperature and UV radiation mostly shapes in the system dynamics.
The recovery time of the system after phage removal can be quantitatively estimated by measuring the relative distance \(\Delta (t)\) between trajectories of the perturbed (i.e. with phage-killing) and the original systems in the model, which we defined as
$$\begin{aligned} \Delta (t)=\sqrt{4\Big ( \frac{S_k(t)-S(t)}{S_k(t)+S(t)}\Big )^2+4\Big ( \frac{I_{1k}(t)-I_1(t)}{I_{1k}(t)+I_1(t)}\Big )^2+4\Big ( \frac{I_{2k}(t)-I_2(t)}{I_{2k}(t)+I_2(t)}\Big )^2+4\Big ( \frac{P_k(t)-P(t)}{P_k(t)+P(t)}\Big )^2}, \end{aligned}$$
(3)
where \(S_k(t), I_{1k}(t),I_{2k}(t),P_k(t)\) are the species densities in the perturbed system, whereas \(S(t), I_1(t), I_2(t),P(t)\) are the densities in the original system, i.e. the one without phage removal. The end of the recovery period can be defined as \(\Delta <\Delta _{cr}\), where \(\Delta _{cr} \ll 1\) is some fixed critical level. Note that along with the introduced metric (3), usage of other similar metrics is possible.
We estimate the impact of the phage removal on the risk of disease acquisition via the following integral quantity G, which we refer to as the integral impact. The integral impact G describes how the exposure of humans to the phage-free bacteria changes after killing free phages, as compared to the original system:
$$\begin{aligned} G=\int _{t_k}^{t_e} (S_k(t)-S(t))dt, \end{aligned}$$
(4)
where \(t_k\) is the time of implementation of the agrochemical (i.e. the date of phage removal), \(t_e\) is the time at which the system recovers after the perturbation, which is determined by the condition \(\Delta (t_e)=\Delta _{cr}\). The introduced integral impact takes into account both the magnitude of the increase in bacterial density after removal of phage and the duration of exposure, i.e. the length of the recovery period.

The integral impact G, characterising the effect of the killing of free phages by agrochemicals to the risk of disease acquisition, is plotted against the total number of hot hours (defined as hours with \(T > 35 \ ^\circ C\)) over the period of the system’s recovery. The non-spatial model is used. Simulations are combined for the eight considered provinces of Thailand. For each province, free phages are assumed to be deactivated (killed) on September 1st each year, within the period 2009–2023. The period of recovery of the system after phage killing, has been estimated using the condition \(\Delta (t)>\Delta _{cr}\). The mathematical expressions for G and \(\Delta (t)\) are given by Eqs. (4) and (3), respectively. We assume \(\Delta _{cr}=0.25\). The unit of G is cells/ml \(\cdot\) day.
Our simulation reveals that G exhibits a pronounced inter-annual variation, for example it is much smaller for the scenario shown in Fig. 7A,C as compared to the situation presented in Fig. 7B,D. We find that the magnitude of G is mainly determined by the temperature, in particular, by the number of hot hours, where \(T > 35 \ ^\circ C\). Figure 8 plots the values of G for different years (2009–2023) and the eight considered Thailand provinces against the number of hot hours during the period of recovery of the system. One can see that for cool seasons (small numbers of hot hours), the effect of killing phages does not increase the exposure to the bacteria, thus implementation of phage-killing agrochemicals at cooler times does not present extra risk of disease acquisition. This is due to the fact that most bacteria exist as lysogens, therefore the absence of free phages will not largely alter the state of the system. On the contrary, for a supercritical number of hot hours, a large proportion of bacteria exist in a phage-free form S. Therefore, setting \(P = 0\) will release them from the control by phage, resulting in a rapid amplification of bacteria S, attaining the carrying capacity level. This signifies that the application of phage-killing agrochemicals at temperatures \(T > 35 \ ^\circ C\) would be potentially risky for agricultural workers.
Role of vertical spatial heterogeneity
In the previous sections, we modelled bacteria-phage dynamics in the surface water of a rice field, using a non-spatial system, assuming a homogeneous environment. Other scenarios of B. pseudomallei-phage regulation are possible, including interaction between microorganisms in soil under a gradient of the temperature, and this, obviously, requires assessing the spatial aspect of the problem. In the current section, we consider B. pseudomallei-phage interactions in soil affected by climate change as well as effects of agricultural practices on the risk of disease acquisition. Note that an important consequence of the presence of a vertical gradient of the temperature is the coexistence of predominately lysogenic and lytic types of infection at different soil layers for each moment of time. Another deviation is that bacteria-phage interactions in soil are not affected by the UV solar radiation, so the major environmental factor, externally forcing the system, is the surface temperature. Finally, we must stress that exhaustive analysis of the spatial model should be done in a separate study, therefore here we only provide a short summary of the main findings.

Simulated annual average densities of susceptible bacteria S, bacteria infected with phage in the lysogenic cycle \(I_1\), bacteria infected with phage in lytic cycle \(I_2\) and the phage P in soil at different depths, as predicted by the spatially explicit model for the Nakhon Phanom province for the period 2009–2023 using historic temperature data. The unit of the densities of bacteria is cell/ml and for phages is pfu/ml.
The outcome of bacteria-phage interaction in soil for the period 2009 to 2023, using historic hourly temperature data, is presented in Fig. 9, constructed for the Nakhon Phanom province. For each year, at each depth, we show the annual average densities of \(S, P, I_1, I_2\). The spatial model predicts the occurrence of two major peaks of the density of susceptible bacteria S: one peak arises near the surface, the second one is located at a depth of approximately 30 cm. The near surface peak exhibits pronounced daily and annual fluctuation, whereas the deeper peak changes on a slower timescale. Similar patterns of vertical profiles of \(S, P, I_1, I_2\), was found in the work of Egilmez and colleagues22, who used a simplified version of the model. The emergence of the deep (the second) peak of S occurs since lysogenised bacteria \(I_1\) experience extra mortality due to lysis during warm periods. In quasi-absence of free phages P at the considered depths (30–40 cm), this gives advantage to the phage-free bacteria S in their competition for resources with \(I_1\). Note that the above-mentioned extra mortality of \(I_1\) due to lysis is generally very small, and it occurs only during relatively short warm periods, where extra heat penetrates the considered depths. This explains the slowness of the emergence of the second peak of S predicted by the model.
Examples of the daily variation of vertical distribution of bacterial and phage densities from our simulation are provided in Supplementary Material SM6. Figure 9 demonstrates that the raise of the surface temperature over the period 2009–2023 results in an overall increase of the abundance of susceptible bacteria in soil. The figure also demonstrates that the development of the deep peak of susceptible bacteria is a slow process, which occurs over a period of 12 years. For the other considered provinces of Thailand, the amplitude and the time of formation of the deep peak of S largely varies from one province to another (we do not show the corresponding figures for brevity).

Predicted annual average vertical distributions for the densities of susceptible bacteria S, bacteria infected with phage in the lysogenic cycle \(I_1\), bacteria infected with phage in the lytic cycle \(I_2\) and the phage P in soil in the space-explicit model for years 2024–2044. The unit of the densities of bacteria is cell/ml and for phages is pfu/ml.
Using our surface temperature forecast from 2024 to 2044, we predicted the dynamics of bacteria and phage density for this indicated period. The corresponding results are shown in Fig. 10 for the Nakhon Phanom. This figure demonstrates that an increase in the surface temperature caused by climate change results in an intensification of both peaks of the density of susceptible bacteria S. The annual average of free phage density in the upper layer of soil follows an increase in the amount of susceptible bacteria. Note that unlike the bacteria-phage dynamics in the surface water, the spatial model does not predict eventual extinction of phage. On the contrary, the total amount of phage towards year 2044 is predicted to increase by approximately a factor of 2. This occurs due to the gradual increase in the amount of susceptible host for the phage because increase in UV radiation would have no effect on either bacterial or phage mortality. The amount of infected bacteria in the lysogenic cycle slightly decreases with time, whereas the numbers of infected bacteria in the lytic cycle show a slight increase. Our simulation of bacterial and phage dynamics in other provinces shows similar patterns of the forecast (we do not show here the corresponding graphs for brevity).
Finally, we modelled the influence of agricultural practices on the risk of disease acquisition in the system with space. Namely, we explored effects of two following agricultural activities: (i) over-turning and mixing the top soil layer due to plowing or shoveling, and (ii) killing the free phages as a result of usage of agrochemicals. In the model, for soil mixing, we assume that each year on May 1st the upper soil layer (from the surface to a depth of 40 cm) becomes mixed. Therefore, on the day of mixing, the densities \(S, P, I_1, I_2\) stay uniform in the mixed layer, with the values corresponding to the spatially averaged values in the considered layer before mixing. We found that soil mixing may have a significant impact on bacterial densities at depths between 20cm and 40 cm. On the other hand, the density of susceptible bacteria near the surface recovers quickly and reaches similar levels to those without mixing. We also found that the density of bacteria decreases at depths of 20 cm to 40 cm due to soil mixing. Soil mixing hinders the development of the deep maximum of susceptible bacteria in soil. Examples of dynamics of the system with and without mixing for the Nakhon Phanom are provided in Supplementary Material (SM7).
To model effects of agrochemicals, we eliminated free phages P in the vertical layer from the surface up to a depth of 5 cm on September 1st each year. Simulation reveals similar outcomes of phage killing as in the homogeneous model, see the previous section. In particular, the application of phage-killing agrochemicals strongly depends on the ambient temperature. In warmer years, it takes approximately two weeks for the phages to regain their regulation of bacteria, while in colder years, it takes about three weeks. On the other hand, the increase of the density of susceptible bacteria due to phage removal in warmer years can be much higher than in colder years. Therefore, application of agrochemicalsis not recommended during periods of high temperatures (\(T > 35 \ ^\circ C\)). Examples of dynamics of the system after removal of free phages in the top 5cm of soil are provided in Supplementary Material (SM7).