The sensitivity and significance analysis of parameters in the model of pH regulation on lactic acid production by Lactobacillus bulgaricus

Background The excessive production of lactic acid by L. bulgaricus during yogurt storage is a phenomenon we are always tried to prevent. The methods used in industry either control the post-acidification inefficiently or kill the probiotics in yogurt. Genetic methods of changing the activity of one enzyme related to lactic acid metabolism make the bacteria short of energy to growth, although they are efficient ways in controlling lactic acid production. Results A model of pH-induced promoter regulation on the production of lactic acid by L. bulgaricus was built. The modelled lactic acid metabolism without pH-induced promoter regulation fitted well with wild type L. bulgaricus (R2LAC = 0.943, R2LA = 0.942). Both the local sensitivity analysis and Sobol sensitivity analysis indicated parameters Tmax, GR, KLR, S, V0, V1 and dLR were sensitive. In order to guide the future biology experiments, three adjustable parameters, KLR, V0 and V1, were chosen for further simulations. V0 had little effect on lactic acid production if the pH-induced promoter could be well induced when pH decreased to its threshold. KLR and V1 both exhibited great influence on the producing of lactic acid. Conclusions The proposed method of introducing a pH-induced promoter to regulate a repressor gene could restrain the synthesis of lactic acid if an appropriate strength of promoter and/or an appropriate strength of ribosome binding sequence (RBS) in lacR gene has been designed.


Background
Lactobacillus delbrueckii subsp. bulgaricus has been widely applied in diary industry, especially as a starter for yogurt production. This microorganism is facultative anaerobic and can produce lactic acid. During the later period of fermentation and storage of yogurt, the production of lactic acid is dominated by L. bulgaricus [1]. The excessive production of lactic acid during storage makes the yogurt taste too sour and this phenomenon is called postacidification.
At industrial level, several attempts have been made to slow down post-acidification, such as using devoid of L. bulgaricus, selecting the weak post-acidification L. bulgaricus, pasteurization after fermentation and so on [2][3][4]. However, these methods either control the postacidification inefficiently or kill the probiotics in yogurt. Using genetic technology to modify the bacteria has become another option to decline post-acidification. Due to lacking molecular tools, mainly caused by the absence of a reliable transformation procedure, our understanding of the physiology and genetics of L. bulgaricus is still limited [5].
In L. bulgaricus, the absorption of lactose is processed by lactose/galactose antiport transport system [6]. Firstly, lactose is transported into cells by lactose permease encoded by lacS gene. Then the lactose is decomposed into glucose and galactose by β-galactosidase encoded by lacZ gene. The free galactose is pumped out of cells or stored in cells in the form of macromolecule i.e. carbohydrate gum. Glucose turns into pyruvate through glycolysis, and the glycolysis is transformed into lactic acid by the catalysis of lactate dehydrogenase [7][8][9]. The lacS gene and the lacZ gene constitute a lacSZ operon. In the downstream of lacZ gene, there is a lacR gene encoding a repressor which makes the lacSZ operon induced by lactose. However, in L. bulgaricus, the lacR gene has lost regulatory function due to the insertion of some gene fragments, resulting in constitutive expression of lacSZ operon [10].
Due to our limited knowledge, most attempts at the genetic level focused on changing the activity of one significant enzyme at a time in the metabolism of lactic acid. Druesne altered the Histidine codon to Alanine codon at the 552th locus of lactose permease, resulting in the reduction of the enzyme's activity [11]. Early in 1990, Molletobtained the mutant L. bulgaricus without βgalactosidase activity by spontaneous deletion [12]. The new aforementioned organisms did produce less lactic acid but they had trouble in growing in milk independently due to the lack of energy. Adams selected two coldsensitive mutants of the β-galactosidase from L. bulgaricus by using the expression system of E. coli [13]. However, their later research turned back to directly mutate and screen L. bulgaricus since it was difficult to transform with L. bulgaricus.
So far, Lactobacillus delbrueckii subsp. bulgaricus ATCC 11842 has been reported successfully transformed by electroporation, although with very low reproducibility and efficiency [5,14]. Several pH-induced promoters from Lactococcus lactis have been demonstrated, such as rcfB [15], P1 and P170 [16]. Combining the two points mentioned above, we came up with a new method to build a pHinduced promoter with a repressor gene controlling the production of lactic acid. Thus, we could turn on the switch at appropriate time. Here we analyse the parameters in the kinetic model to investigate the regulation effect of pH-induced promoter on lactic acid.

The reduced metabolic network
Here we present a model of the production of lactic acid, which includes the important enzymes, lactose permease, β-galactosidase, lactate dehydrogenase and related regulation genes lacS, lacZ, lacR, as shown in Figure 1. The process of glycolysis is reflected in one reaction, adopting the assumption that there is a constant flux from glucose to pyruvate. Pyruvate is partly catalysed to lactic acid, and the other two products, acetyl-CoA and butanedione, are not taken into account.

The mathematical model
The mathematical model describes how the extracellular pH value influences the production of lactic acid. The ordinary differential equations in the model are based on Michaelis-Menten equations [18] and fundamental kinetic principles. The variables are described in Table 1. The model assumptions are as follows.
(1) Each entity, described in Table 1 is given as total amount in the population.
(2) The model does not take energy metabolism into consideration.
(3) All reactions are modelled by mass action principles, except for enzyme reactions which obey Michaelis-Menten equation, and transcription which obeys saturation kinetics. (4) The level of translation product is assumed to be proportional to mRNA levels according to Dockery and Keener [19]. Figure 1 A schematic of metabolic and gene regulation network model of the production of lactic acid. The LacS, LacZ and LDH represent lactose permease, β-galactosidase, and lactate dehydrogenase respectively. When the extracellular pH reduces to 5.5, the pH-induced promoter (we chose rcfB promoter as an example) will trigger the express of LacR. The accumulation of LacR will combine with the operator sequence upstream of lacS gene, inhibiting the expression of LacS and LacZ [17].

The metabolism of lactose
As mentioned before, the metabolism of lactose is described by Michaelis-Menten equation. Since lactose transport is reversible, a term was included to count for lactose efflux dependent on the internal lactose concentration [20]. The V max for lactose permease and β-galactosidase is associated with the concentration of enzymes, so we use K cat and the concentration of enzymes to estimate the V max . Since the lac operon is usually induced by lactose, we assume that the repressor LacR protein could react with lactose. Thus we need to take the reduction into consideration when evaluating the concentration of intracellular lactose and LacR protein. As for the utilization of sugar, L. bulgaricus belongs to the type of homofermentation. Theoretically, this microorganism could totally convert glucose to lactic acid, but the conversion of glucose is about 80%-90% in practice [21]. Therefore, the estimated proportion of pyruvate transformed into lactic acid is 80%.

The gene regulation
The transcription and decomposition of lactose are tightly regulated at the genetic level: production of the enzymes can be decreased at the transcriptional level by regulatory protein LacR binding at appropriate DNA sites. We set M to stand for the transcript concentration of lacR gene. Since lacS and lacZ are located in the same operon, it is assumed that both genes are transcribed at same rate. We used N to describe the concentration of the two genes' mRNA. When describing the equations, the mRNAs are assumed to be transcribed at two distinct rates: basal one when there is no regulation and a higher/lower rate when being regulated [22].

The transcription of lacR gene
For M, it is described as function (1). The basal rate is given by V 0 and increased at a constant rate V 1 which is regulated by the pH value. GR stands for the repressor gene concentration. The degradation of the mRNA occurs at the rate of d M . It is also diluted due to cell growth rate at μ. We assumed the decrease of mRNA is first-order proportion to the concentration of M.
As the mechanisms of pH-induced promoters have not been fully elucidated, we included a generic pH-dependent switch F(pH) which turns on when pH value is below the threshold. The switch takes the form of the smoothed step function [23], where 5.5 represents the threshold pH level around which the switch occurs [15]. n indicates the steepness of the smooth switch function.
Since lactic acid is the main product during fermentation, we assumed pH value is all controlled by the concentration of lactic acid. Besides, we needed to take the buffer capacity of milk medium into account. Furthermore, we assumed that all the lactic acid could spread to the medium even when the extracellular concentration was pretty high. We tested the pH value of milk medium when different amount of lactic acid was added. Three parallel experiments were taken and we fitted a function as the form of function (3). To the pH data, c 0 , c 1 and c 2 are constants. This function was then fed into the model via the switch function F(pH). 2 (3)

Transcription of lacS and lacZ gene
For N, it is described as function (4). The basal rate is given by V 2 and the transcription is inhibited by LacR. For the part to simulate the regulation rate of LacR, T max stands for the maximal rate of lacS and lacZ gene transcription, S represents the sensitivity of lacS and lacZ gene transcription to lactose permease and β-galactosidase [24]. G is the concentration of transcription gene. Again, we assumed first-order degradation rate at d N and dilution rate at μ of the mRNA.

Cell growth
In milk medium, lactose is the main carbon source, so it is assumed that the cell growth rate is dependent on the extracellular lactose concentration. The product, lactic acid, which declines the pH value of milk medium performs inhibition on cell growth. The function to describe the production inhibition is the same as Concepcion presented [21]. In function (5), K s is the Monod constant for growth in extracellular lactose, and μ max is the maximum specific growth rate. K LA is the maximum  initial lactic acid concentration in which the microorganism growth is completely inhibited.
Combining the metabolism of lactose with gene regulation All the definitions of the parameters are shown in Table 2.
Combining with the enzyme reactions, gene regulations and cell growth discussed above, the process described in Figure 1 could be represented by: The final metabolic model is given by: We set stoichiometric constant of R 4 with two since two molecules of pyruvate are formed from one molecule of glucose. For the stoichiometric constant of R 5 , it stands for the proportion of pyruvate which is converted to lactic acid.

Data estimated
Firstly, we did not take lacR gene into consideration. This is just the case of wild type L. bulgaricus which lacSZ operon is constitutive expression. Since most of the parameters in this reduced model have been reported, V 2 is the only one needs to be modified. We used SBToolbox in Matlab to set up the reduced model and simulated the metabolism process [30]. We compared the concentration changes of extracellular lactose and lactic acid within 8 hours with the data from Fatama, where the conditions used were 43°C, 8% of skim milk concentration (without pH value control) and 4% of inoculum ratio [31]. After that, we added lacR gene. This time we did not have sufficient data for full parameterization of the model. However, we did have enough information to estimate their relative sizes and give the qualitative nature of our investigations. The parameter values are listed in Table 2.

Sensitivity and significance analysis of parameters
We still used SBToolbox to perform local sensitivity analysis and Sobol's method for global sensitivity analysis. The local sensitivity analysis investigates how small changes in a single parameter value could affect the model output. The method is based on the partial differentiation of the output with respect to the input parameters [32]. Herein, the partial differentiation is evaluated numerically by introducing a 1% increment from the specific parameter value. We chose Sobol sensitivity analysis method to calculate global sensitivity since the sensitivity index quantified the overall effects of a parameter, in combination with any other parameters, on the model output [33]. The number of simulation to carry out is 10000 times. We first set the number as 1000 according to Schmidt [30]. Since the results are varied among different simulations, we increased the number to 10000 and got stable outcomes. During the above sensitivity analysis, we chose lactic acid concentration as the model output.
After figuring out the sensitivity parameters to the model output, we carried out different simulations with each sensitivity parameter at two levels, a relative low value and a relative high value. The metabolism of lactic acid and changing of pH values was carefully considered. Those results indicate how pH values affect lactic acid production under different situations and help us to understand how to design LacR regulation in biology experiments in future.

Reduced model without lacR gene
The results of modelled concentration changing of extracellular lactose and lactic acid are shown in Figure 2, along with the experimental data. The model fits the experimental data well. We have obtained the value of V 2 , which indicates the basal rate of lacSZ gene transcription. The relative size of fitting value (V 2 = 3.5 h -1 ) is determined according to Gustafsson's [24] report.

Sensitivity analysis of parameters
The results of local sensitivity analysis and Sobol sensitivity analysis are shown in Figure 3. The parameters with sensitivity from high to low in local sensitivity analysis are followed as T max , GR, K LR , S, V 2 , d M , V 1 , u max , d LR  and V 0 . The other parameters perform little effect on the production of lactic acid. The reason is that when one of those parameters increased, it would only affect the production rate of the output but have no effect on the total amount of lactic acid. Thus, they are not sensitive parameters in this method. As for Sobol sensitivity analysis, the parameters with sensitivity from high to low is T max , K m·LP , K LZ , K LR , GR, K m·LYC , V 0 , S, d LS , V 1 , d N , V 2 and d LR . The parameter K m·LP shows the second most sensitivity in Sobol's method, just as Druesne has proven, the reduction of lactose permease's activity would greatly decline the production of lactic acid [11].

Significant analysis of parameters
In sensitivity analysis, sensitivity does not mean importance, since sensitive parameters are always fixed inherently and unadjustable. We need to combine the sensitivity results with operability, especially for the design of biology experiments. Although T max , K LZ , V 2 , d LS , d LR , GR, d M , d N and u max perform high sensitivity in both methods, these parameters are usually inherently. K m·LP value is adjustable, but to decrease its value would generate fatal disadvantage on the growth of the bacteria. So we do not include this parameter into consideration. Due to the same reason, we exclude K m·GLYC , which is the Michaelis constant for glycolysis. As for the parameter S, which means sensitivity of lacSZ gene transcription to LacS and LacZ proteins, it is difficult to measure and adjust in experiments. Therefore, the parameters which exhibit great influence on the production of lactic acid in both methods are K LR , V 0 and V 1 , with K LR representing RBS strength of lacR gene, and V 0 together with V 1 reflecting characters of the pH-induced promoter.
For the three parameters, the relative high levels of K LR , V 0 and V 1 are 1128 h -1 , 1 h -1 and 21 h -1 respectively. The relative low levels of K LR , V 0 and V 1 are 564 h -1 , 0.1 h -1 , 10.5 h -1 respectively. The results are shown in Figure 4. From picture A, we could figure out that with the increase of promoter strength, the production rate and amount of lactic acid would be greatly decreased. However, for the relative high value of V 1 , the finial pH value only declines to 4.83 which is higher than casein isoelectric point. Thus we need to find an appropriate promoter strength between the relative high and low values when it is turned on by pH value. From picture B, we could draw the conclusion that if the pHinduced promoter could be efficiently turned on when pH decreases to its threshold, leaking expression of the promoter would have little effect on the production of lactic acid. In picture C, the increase of K LR leads to strong decline of both the production rate and total amount of lactic acid. By comparing the results of relative high level K LR with the relative high level V 1 in picture A, the two parameters almost perform the same Figure 2 Simulation of batch fermentation conversion of lactose to lactic acid using wild type L. bulgaricus. The initial conditions were 43°C 8% of skim milk concentration (without pH value control) and 4% of inoculum ratio. The dots and triangles indicated the measured concentration of lactose and lactic acid in the medium, respectively, every two hours. The curves were the simulated results of lactose and lactic acid concentration in 8 hours.  influence on the production of lactic acid. This is due to the reason that the two parameters both affect the concentration of LacR. Therefore, we could adjust the strength of pH-induced promoter and/or RBS strength of lacR gene to make our switch perform well in controlling lactic acid production.
Other two important features for the promoter are the threshold pH level around which it would start and whether the promoter would keep turning on when the pH value decreased to 3.5 or lower. Akyol [15] and Madsen [16] only described the status of pH-induced promoters along with the pH declined to 5.5. Further experiments are needed to test the characteristics of the promoters. This model neglect the effect of lower pH on the promoter, assuming that it would still turn on with the decrease of pH value.

Conclusions
In conclusion, we propose a new method to control the production of lactic acid by L. bulgaricus which is to build a pH-induced switch (promoter) to turn on repressor gene. The proposed method overcomes the disadvantage in bacterial growth by directly changing one of the enzyme related to lactic acid metabolism. Then, we build a reduced model of pH-induced promoter regulation on the production of lactic acid and adjust the model with data from wild type L. bulgaricus. After that we carry out sensitivity analysis of the parameters and figure out three significant ones in the model. To make our switch work well, we need to find an appropriate strength of promoter and/ or an appropriate strength of RBS in lacR gene. The values of those two parameters should between the high level and low level we have set in the analysis, so that when the pH value declines to the promoter's starting point, the lacR gene could express moderately to inhibit the production of lactic acid.