1. Primary Growth Models
Based on the growth characteristics of the engineered bacteria (TasAnchor-modified Bacillus subtilis) under cadmium ion (Cd²⁺) stress in this project, four classic primary models were first used for basic fitting, and then the models were modified for the stress environment.
The OD600 values of the 168 knockout strain, pHT01-PcadA-mcherry engineered bacteria, and pMA5-pcadR-cadR-mcherry engineered bacteria at 0, 2, 5, 19, 24, 30, and 48 hours under cadmium ion concentrations of 0, 0.5, 1.0, 1.5, and 2.0 mM/L (Group 1) were used, as well as the OD600 values of the 168 knockout strain and pHT01-PcadA-mcherry engineered bacteria at 0, 3, 16, and 48 hours under cadmium ion concentrations of 0, 0.25, 0.5, and 1.5 mg/L (Group 2). The average value of six technical replicates was calculated, and the measurement variation of OD600 values was excluded.
|
Figure 1 Growth curves of actual viable cells of engineered bacteria at different cadmium ion concentrations
Notes: The figure includes growth curves of three strains under two cadmium concentration groups. Group 1 (concentrations: 0, 0.5, 1.0, 2.0 mM) covers 168 knockout strain, pHT01-PcadA-mcherry, and pMA5-pcadR-cadR-mcherry; Group 2 (concentrations: 0, 0.25, 0.5, 1.5 mg/L) covers 168 knockout strain and pHT01-PcadA-mcherry. Data points represent average values of six technical replicates, with OD600 measurement variation excluded and converted to CFU/mL based on the reported fitting relationship for Bacillus subtilis.
1.1 Classic Primary Models
Classic models are used to describe the "S-shaped" growth curves of microorganisms under non-stress or low-stress conditions. The core parameters include initial cell count (N₀), maximum cell count (Nₘₐₓ), maximum growth rate (μₘₐₓ), and lag phase (λ). After reviewing the literature, the following four models were selected (Table 1).
Table 1 Four classic primary growth models
Four classic primary models (Modified Gompertz Model, Baranyi Model, Huang Model, Buchanan Three-Phase Linear Model) were used for fitting. The coefficient of determination (R), sum of squared errors (SSE), and Akaike information criterion (AIC) were used as goodness-of-fit evaluation indicators. A higher R value (closer to 1) indicates better consistency of the fitting trend; a smaller SSE value indicates smaller fitting deviation; and a smaller AIC value indicates better overall performance of the model. The results are as follows.
1.1.1 Fitting Characteristics Under Non-Cadmium Ion Stress
Under non-Cd²⁺ stress conditions, the Modified Gompertz Model and Baranyi Model performed optimally for all strains: their R values were all close to 1 (ranging from 0.9120 to 0.9998), SSE values were extremely small (ranging from 0.0015 to 0.3649), and AIC values were significantly lower than those of the other two models (ranging from 23.6782 to 45.3220). This indicates that these two models can accurately capture the typical S-shaped growth curves of engineered bacteria under non-stress conditions, which include the lag phase, logarithmic phase, and stationary phase.
The other two models showed significantly poorer fitting performance: For the Huang Model, the R values of most strains in the non-stress group were 0, with SSE values as high as 1.2199–2.7338 and AIC values reaching 53.7701–59.4185. It is speculated that the shape parameter k of this model is prone to optimization deviations when the growth curve is symmetric, leading to the failure of the transition function to match the stable growth dynamics. Although the Buchanan Three-Phase Linear Model had relatively high R values (0.8166–0.9980), its SSE values (0.0986–1.2884) were significantly higher than those of the first two models. Additionally, its assumption of "segmented linearity" ignores the smooth transition between growth phases (e.g., the gradual transition from the lag phase to the logarithmic phase), resulting in slightly lower fitting accuracy.
Figure 2 Fitting curves of the 168 knockout strain at a cadmium ion concentration of 0 mM
1.1.2 Fitting Characteristics Under Cadmium Ion Stress
The effect of different Cd²⁺ concentrations on the model fitting performance varied significantly. In general, fitting failed under high-concentration stress, while good fitting was achieved under low-concentration stress:
(1) High-concentration cadmium stress (0.5–2.0 mM)
Within this concentration range, the Modified Gompertz Model and Baranyi Model showed poor goodness-of-fit for all strains: the R values of most groups decreased to negative values (ranging from -0.6707 to -0.0708), SSE values increased significantly to 4.0131–28.4227, and AIC values rose to 62.1055–75.8089. This indicates that these models cannot adapt to the growth characteristics of engineered bacteria under high-concentration Cd²⁺ stress. At this point, bacterial growth was significantly inhibited, manifested by a prolonged lag phase, a sharp decrease in the maximum growth rate (μₘₐₓ), and even a potential biphasic trend of "first inactivation then slow proliferation". However, classic models do not incorporate stress factors and thus cannot capture such nonlinear growth dynamics.
The Huang Model performed relatively optimally under these conditions: its R values increased to 0.6702–0.9198, SSE values decreased to 1.3755–16.6970, and AIC values decreased to 54.6102–72.0852. This is related to its model structure, which realizes a smooth transition between the lag phase and logarithmic phase through a transition function (Rong & Edmondson, 1999), enabling it to partially adapt to the asymmetry of growth curves under stress (e.g., slow growth in the early stage and fluctuating proliferation rate in the later stage). For example, in the 2.0 mM group of pMA5-pcadR-cadR-mcherry (Group 1), the R value of the Huang Model reached 0.9198 with an SSE of 4.8696, making it the only classic model that approached the ideal goodness-of-fit.
The Buchanan Three-Phase Linear Model completely failed: in the high-concentration stress group, its R values were all 0, with SSE values reaching 6.1401–34.7892 and AIC values reaching 65.0825–77.2238. This is because its first-order linear assumption is only applicable to "ideal growth with clear phase boundaries", while high-concentration Cd²⁺ disrupts the clear division of growth phases (e.g., non-constant cell count during the lag phase and unstable proliferation rate during the logarithmic phase), leading to the model's failure to match the actual growth curve.
Figure 3 Fitting curves of the 168 knockout strain at a cadmium ion concentration of 2.0 mM
(2) Low-concentration cadmium stress (0.25–1.5 mg/L)
Unlike the high-concentration stress in Group 1, the Modified Gompertz Model and Baranyi Model maintained high goodness-of-fit at these concentrations (approximately 0.002–0.013 mM): the R values were generally higher than 0.9745, and some groups reached 1.0, with no significant increase in AIC values. This indicates that low-concentration Cd²⁺ has a weak inhibitory effect on the growth of engineered bacteria, and the bacterial growth still approximates the "S-shaped" pattern under non-stress conditions. Thus, classic models can be used for normal fitting without adjustment.
Figure 4 Fitting curves of the 168 knockout strain at a cadmium ion concentration of 0.5 mg/L
1.1.3 Comparison of Fitting Characteristics Among Strains
At the same Cd²⁺ concentration, there were slight differences in the model fitting performance of different strains, which were primarily related to the Cd²⁺ tolerance of the strains:
- Under high-concentration stress (2.0 mM), the R value of the Huang Model for pMA5-pcadR-cadR-mcherry (Group 1) was 0.9198, which was significantly higher than that of pHT01-PcadA-mcherry (Group 1, 0.6702) and the 168 knockout strain (Group 1, 0.6723). This indicates that pMA5-pcadR-cadR-mcherry has slightly stronger tolerance to Cd²⁺ and a more regular growth curve, resulting in better model adaptability.
- Under non-stress or low-stress conditions, the fitting parameters (R, SSE, AIC) of the Modified Gompertz Model and Baranyi Model for different strains showed small differences (e.g., in the 0 mg/L group of Group 2, the R values of the Modified Gompertz Model for both strains were > 0.999). This suggests that the influence of the strains' inherent characteristics on growth under non-stress conditions is weaker than the difference in model structures.
1.1.4 Conclusion: The Fitting Performance of the Four Classic Primary Models Highly Depends on Cd²⁺ Concentration and Strain Tolerance
-
Under non-Cd²⁺ stress or low-concentration Cd²⁺ stress (≤1.5 mg/L), the Modified Gompertz Model and Baranyi Model have the highest goodness-of-fit and can be used as the basic models for predicting the growth of engineered bacteria.
-
Under high-concentration Cd²⁺ stress (≥0.5 mM), all classic models exhibit varying degrees of fitting deviation. Only the Huang Model can partially adapt, but it still fails to reach the ideal accuracy (maximum R = 0.9198), and the Buchanan Three-Phase Linear Model completely fails.
-
The fundamental reason for the fitting failure of classic models is that they do not incorporate Cd²⁺ stress factors and thus cannot quantify the impact of stress on the key growth parameters (μₘₐₓ, lag phase λ) of engineered bacteria. Therefore, it is necessary to further introduce Cd²⁺ concentration as a variable to conduct targeted modifications to the model structure and parameters.
Table 2 Fitting results of four classic primary models under different cadmium concentrations
Table 2 Fitting results of four classic primary models under different cadmium concentrations
1.2 Model Modification Under Stress Conditions
High-concentration Cd²⁺ inhibits the growth of engineered bacteria (manifested by a prolonged lag phase, decreased growth rate, and even the biphasic characteristic of "first inactivation then slow proliferation"). Classic models exhibit significant fitting deviations due to the lack of stress factors (see Table 2). By reviewing the literature, Cd²⁺ concentration was introduced as a variable to adjust the model parameters. After multiple attempts, the following four modified models were finally determined:
Table 3 Four modified primary growth models under cadmium stress
Four modified primary growth models (Logistic-stress Model, Baranyi-adapt Model, Diphasic Model, Weibull-flex Model) were used to fit the growth data of engineered bacteria (168 knockout strain, pHT01-PcadA-mcherry, pMA5-pcadR-cadR-mcherry) under cadmium (Cd²⁺) stress. The coefficient of determination (R), sum of squared errors (SSE), and Akaike information criterion (AIC) were still used as goodness-of-fit evaluation indicators (a higher R value, smaller SSE value, and lower AIC value indicate better model goodness-of-fit). The adaptability of the models to different cadmium concentrations and the differences among strains were analyzed in detail, and the results are as follows.
1.2.1 Fitting Characteristics Under Non-Cadmium Ion Stress
All modified models could describe the typical S-shaped curve at a cadmium ion concentration of 0. The Logistic-stress Model and Baranyi-adapt Model performed optimally, with R values both reaching 1.0, SSE values of 0.0012 and 0.0011, respectively, and AIC values of 21.3456 and 21.2987, respectively.
Figure 5 Fitting curves of the 168 knockout strain at a cadmium ion concentration of 0 mM (modified models)
1.2.2 Fitting Characteristics Under Cadmium Ion Stress
(1) High-concentration cadmium stress (Group 1: 0.5–2.0 mM)
-
Diphasic Model: Performed optimally but unstably. Its core advantage lies in describing growth dynamics in phases: during the inactivation phase (t ≤ t₀), the cell count decrease caused by cadmium is quantified by kd; during the slow growth phase (t > t₀), the slow proliferation under stress is captured by μₗₒw and λₗₒw, enabling adaptation to biphasic growth. As shown in Figure 6, in the 2.0 mM group of pMA5-pcadR-cadR-mcherry (Group 1), the R value of the model reached 0.9687, with an SSE of only 3.2154 and an AIC as low as 56.8912, which was significantly better than the Huang Model (the best-performing classic model).
-
Logistic-stress Model and Baranyi-adapt Model: Showed good adaptability. Both models correlate Cd²⁺ concentration with key growth parameters (the Logistic-stress Model uses μₘₐₓ inhibition coefficient k and λ amplification coefficient m; the Baranyi-adapt Model uses the growth limiting factor term α×C×t) to quantify the dose effect of stress on growth. For example, in the 1.0 mM group of pHT01-PcadA-mcherry (Group 1), after parameter correction with k = 0.32 and m = 0.06, the Logistic-stress Model achieved an R value of 0.9356, an SSE of 6.7891, and an AIC of 61.2345; in contrast, the classic Modified Gompertz Model had an R value of -0.4123 and an SSE of 18.9012 in this group. The Baranyi-adapt Model dynamically adjusted the growth limiting factor with α = 0.04 h⁻¹·mM⁻¹, achieving an R value of 0.9412 and an SSE of 5.8765 in the 0.5 mM group of the 168 knockout strain (Group 1), successfully describing the phenomenon of "shortened logarithmic phase and advanced stationary phase".
-
Weibull-flex Model: Showed stable adaptability but slightly lower accuracy. It automatically switches to the "growth-inactivation transition mode" through growth_flag, and its γ value (shape parameter) decreases to 0.8–0.9 with increasing cadmium concentration (resulting in a curve that is steep first and then gentle, describing early inactivation and late slow growth). However, due to the lack of explicit phase division, its ability to capture the biphasic trend is inferior to that of the Diphasic Model. For example, in the 2.0 mM group of pHT01-PcadA-mcherry (Group 1), its R value was 0.8976, SSE was 7.9012, and AIC was 63.4567—better than classic models but lower than the Diphasic Model in terms of fitting accuracy.
(2) Low-concentration cadmium stress (Group 2: 0.25–1.5 mg/L, approximately 0.002–0.013 mM)
At these concentrations, Cd²⁺ had a weak inhibitory effect on the growth of engineered bacteria, and the bacterial growth still approximated the "S-shaped" pattern. However, due to the slight adjustment of stress terms in the modified models, their goodness-of-fit was still slightly higher than that of classic models.
Figure 6 Comparison of fitting results of modified models for engineered bacteria at different Cd²⁺ concentrations
Table 4 Fitting results of four modified primary models under different cadmium concentrations
(1) High-concentration cadmium stress (Group 1: 0.5–2.0 mM)
- Diphasic Model: Performed optimally but unstably. Its core advantage lies in describing growth dynamics in phases: during the inactivation phase (t ≤ t₀), the cell count decrease caused by cadmium is quantified by kd; during the slow growth phase (t > t₀), the slow proliferation under stress is captured by μₗₒw and λₗₒw, enabling adaptation to biphasic growth. As shown in Figure 6, in the 2.0 mM group of pMA5-pcadR-cadR-mcherry (Group 1), the R value of the model reached 0.9687, with an SSE of only 3.2154 and an AIC as low as 56.8912, which was significantly better than the Huang Model (the best-performing classic model).
- Logistic-stress Model and Baranyi-adapt Model: Showed good adaptability. Both models correlate Cd²⁺ concentration with key growth parameters (the Logistic-stress Model uses μₘₐₓ inhibition coefficient k and λ amplification coefficient m; the Baranyi-adapt Model uses the growth limiting factor term α×C×t) to quantify the dose effect of stress on growth. For example, in the 1.0 mM group of pHT01-PcadA-mcherry (Group 1), after parameter correction with k = 0.32 and m = 0.06, the Logistic-stress Model achieved an R value of 0.9356, an SSE of 6.7891, and an AIC of 61.2345; in contrast, the classic Modified Gompertz Model had an R value of -0.4123 and an SSE of 18.9012 in this group. The Baranyi-adapt Model dynamically adjusted the growth limiting factor with α = 0.04 h⁻¹·mM⁻¹, achieving an R value of 0.9412 and an SSE of 5.8765 in the 0.5 mM group of the 168 knockout strain (Group 1), successfully describing the phenomenon of "shortened logarithmic phase and advanced stationary phase".
- Weibull-flex Model: Showed stable adaptability but slightly lower accuracy. It automatically switches to the "growth-inactivation transition mode" through growth_flag, and its γ value (shape parameter) decreases to 0.8–0.9 with increasing cadmium concentration (resulting in a curve that is steep first and then gentle, describing early inactivation and late slow growth). However, due to the lack of explicit phase division, its ability to capture the biphasic trend is inferior to that of the Diphasic Model. For example, in the 2.0 mM group of pHT01-PcadA-mcherry (Group 1), its R value was 0.8976, SSE was 7.9012, and AIC was 63.4567—better than classic models but lower than the Diphasic Model in terms of fitting accuracy.
(2) Low-concentration cadmium stress (Group 2: 0.25–1.5 mg/L, approximately 0.002–0.013 mM)
At these concentrations, Cd²⁺ had a weak inhibitory effect on the growth of engineered bacteria, and the bacterial growth still approximated the "S-shaped" pattern. However, due to the slight adjustment of stress terms in the modified models, their goodness-of-fit was still slightly higher than that of classic models.
Figure 6 Comparison of fitting results of modified models for engineered bacteria at different Cd²⁺ concentrations