The model part provides us an opportunity to attempt some interdisciplinary methods and one of the most important features of modeling is that we can design experiments reasonably and have a better knowledge of some quantitative features of the biological system after modeling. Interestingly, models can also help us to discover something novel and achieve goals that can’t be fulfilled with only traditional experiments.

In our project, the engineered bacteria (E. coli modified by design) limit the growth of the Propionibacterium acnes (P. acnes) in two ways: the engineered bacteria can secrete bacteriocins, which directly kill the Propionibacterium acnes; at the same time, the engineered bacteria compete with the Propionibacterium acnes for nutrient fatty acids.

We've developed two models to explain how our co-cultural system functions and how we accelerate fatty acids consumption of engineered bacteria. These are the growth simulator model and the flux accelerator model which consists of two parts.


1. The growth simulator model predicts the trend of the change of our coculture system and gives feedback to how to modify our product, considering it as a whole.
2. The first part of The flux accelerator model reflects the metabolism of E. coli under different oxygen conditions to a certain extent.
3. The second part of The flux accelerator simulates a more complex metabolic network and successfully selects enzymes which have the most potential to promote the decomposition of higher fatty acids.

Fig. 1 Model overview

In the growth simulator model, we not only explain how bacteria relieve acne on the skin surface, but also predict changes of the number of Propionibacterium acnes, which helps us to consider more influential factors that have not been considered in wet experiments. At the same time, this model simulates the final design of our project, reflecting it as a whole, which theoretically complements the experiments.

In the flux accelerator model, we first build an FBA model based on fatty acid β-Oxidation to analyze the enzymes that we need to overexpress to improve the fatty acids consumption ability of the engineered bacteria. Then, we use an existing molding tool to build a more complex FBA model for a more satisfying result. Deciding which enzymes had most potential to promote the decomposition of fatty acids saves us a considerable amount of work on the less effective enzymes.

2.The growth simulator model


In the growth simulator model, we include Propionibacterium acnes and engineered bacteria into the same system, and construct an ordinary differential equation model simultaneously considering the processes of fatty acid consumption and bacteriocin secretion. The model reasonably explains the mechanism of the two pathways of growth inhibition of Propionibacterium acnes by engineered bacteria. Due to the limit of experimental conditions, the wet experimental group hasn't conducted coculture experiments of the strains. Therefore, the work of our mathematical model group is an extension of the wet experimental work, and at the same time illustrates the feasibility of the acne treatment idea to a great extent.


To simplify our model and emphasize the characters we want to investigate, we bring up following assumptions:

1. The promoter efficiency of bacteriocins is proportional to the concentration of fatty acids in the environment[1], and the rate of the consumption of bacteriocins due to the exposure to the outside environment is considered as a constant.
2. There only exists the competitive effect between the engineered bacteria and the Propionibacterium acnes on fatty acids as well as the killing effect through bacteriocins.
3. To emphasize the effect of fatty acids on engineered bacteria and Propionibacterium acnes, we set the environmental capacity of both as a linear function of fatty acids.
4. The killing effect of bacteriocins on P. acnes remains stable during the whole process.
5. To eliminate the effect of units, we use the percentage of relative initial time as data to fit the parameters.

2.3.Symbols and definitions

Symbols Definitions
S1 Relative number of engineered bacteria
S2 Relative number of Propionibacterium acnes
C1 Concentration of bacteriocins
C2 Concentration of fatty acid

2.4.Parameters and description

Symbols Definitions
r1 Growth coefficient of engineered bacteria
r2 Growth coefficient of P. acnes
p1 Bacteriocin production rate affected by engineered bacteria
p2 Bacteriocin production rate affected by fatty acids
p3 Decrease of bacteriocin per unit time
p4 Increase of fatty acids per unit time
K Ratio of reduction of Propionibacterium acnes to the number of bacteriocins
K1 Rate of fatty acid consumption of engineered bacteria
K2 Rate of fatty acid consumption of Propionibacterium acnes
a1 Growth coefficient 1 of engineered bacteria affected by fatty acids
b1 Growth coefficient 2 of engineered bacteria affected by fatty acids
a2 Growth coefficient 1 of Propionibacterium acnes affected by fatty acids
b2 Growth coefficient 2 of Propionibacterium acnes affected by fatty acids

2.5.Model construction and explanation

We use the following symbol system[2] to represent the interrelationship of the four elements of the system (the bacteria, the Propionibacterium acnes, the bacteriocins and the fatty acids). Considering the indirect effect of fatty acids on bacteriocins, we use dashed lines to indicate indirect facilitation.

Fig. 2 System symbols used in our modeling

We use this symbol system to make a graphical representation of the relationship among engineered bacteria, P. acnes, fatty acids, and bacteriocins as follows.

Fig. 3 Relationship among four elements in our system

This graph shows that both engineered bacteria and P. acnes consume fatty acid, engineered bacteria produce PctA( a kind of bacteriocins), and PctA inhibits P. acnes.

According to the graph, we transform the relationships into a system of ordinary differential equations as a simulation. After collecting data from online databases, we can work out the parameters, then predict the changes of the four parts.

Our ideas for designing each ordinary differential equation are as follows:

1. The first equation: In the modeling process, we have assumed that the growth of the engineered bacteria is limited by Propionibacterium acnes and should roughly show an 'S- shaped growth curve', so we refer to the famous logistic population size change theory to support the growth model of the engineered bacteria.
2. The second equation: We assume that bacteriocins are secreted only by engineered bacteria, so the p1S1 in the equation represents this relationship, and because of the positive effect of fatty acid on promoter efficiency of bacteriocins, we use p2C2 to represent this relationship. Besides, the p3 indicates that the consumption of bacteriocins due to the environment is stable through the whole process.
3. The third equation: The change in the number of Propionibacterium acnes is the same as that of the engineered bacteria in accordance with the Logistic growth model, and there is a corresponding decrease KC1 due to the inhibition of bacteriocins, which for simplicity we assume to be proportional to the concentration of bacteriocins.
4. The fourth equation: We consider that fatty acids are only secreted by human skin at a constant rate p4, while fatty acid consumption is caused by both engineered bacteria and Propionibacterium acnes, and the rate of consumption is proportional to the number of both, and can be described by k1S1 and k2S2 respectively.


Fig. 4 Simulated trend of four elements in our system

Based on the reasonable assumptions, we first gave a set of parameters, whose result can fit with our project well. Bringing them into the Growth Simulator, the final overall effect is shown as follows. And due to the lack of data from real experiments, the horizontal axis representing time doesn’t have real meanings but just indicates the relative time of the change.

By observing the trend of the four main elements in the system over time, we can initially derive the effect of each parameter on each of the four elements.

S1 (engineered bacteria) increases steadily with time, while we can also find a gradual decrease in the rate of growth (which is represented by the slope of the tangent line of the curve) with time. This indicates that in the co-culture system, the number of engineered bacteria increases over time, but this increase is limited by the size of its own population, which is also consistent with the Logistic model assumption.

The number of S2 (P. acnes) climbs with time at the very beginning but then drops suddenly and quickly, which shows the effect of bacteriocins and indicates that engineered bacteria limits the increase of P. acnes by competing for the nutrient substrate fatty acids.

C1 (the concentration of bacteriocins) gradually increases with time, with a strong correlation between its growth and the engineered bacteria. This indicates that the bacteriocins are secreted into the system by the engineered bacteria, and since the rate of secretion is assumed to be constant, a strong correlation ought to exist between the S1 curve and the C1 curve.

The decrease in C2 (concentration of the fatty acid) over time is easily explained by the fact that engineered bacteria and Propionibacterium acnes use fatty acids as nutrients.

2.7.Sensitivity Analysis

In this section, we change the parameters in our model separately by 5%, 10% and 20% with other parameters staying the same, to find the best character to improve our engineered bacteria.

Aimed at relieving acnes, we only choose to analyze the parameters relevant to the change of fatty acids and P. acnes. According to the description of all the parameters, we select r1, p1, p3, K and K1. All the parameters are increased by 5%, 10% and 20% in the following analysis except for p3, which is decreased by different percentages.

In the following figures, the vertical axis represents the value of y, which is calculated with the following definition:

y = y2 - y1

y1 represents the content of material when the content of P. acnes decreases to close to 0 after the change of a parameter. So y1 is a minimal positive value.

y2 represents the original content of material without the change of Parameters at the same time of y1 determined.

Fig. 5 Sensitivity analysis of five selected parameters by 5%

Fig. 5 above shows the difference in the change of percentage of fatty acids and P. acnes before and after the change of parameters shown in the horizontal axis by 5%. We can find that:

1. The changes of r1, p1and K have no significant effect on the relative content of P. acnes. But when p3 or K1 is changed by 5%, the relative content of P. acnes significantly changes, which means that if we manage to enhance the decrease of fatty acids per unit time or the rate of fatty acids consumption by engineered bacteria, we can effectively reduce the number of P. acnes.
2. Except for K1, changes of other parameters have little effect on the relative content of fatty acids. The changes of r1, p1and K even have negative effects because of the increase of the relative number of fatty acids.

Fig. 6 Sensitivity analysis of five selected parameters by 10% and 20%

When it comes to changes by 10%, all parameters have more significant effects on accelerating the reduction of P. acnes content and this kind of trend will continue when parameters are changed by 20%. The effects of changing p3 and K1 are still the most significant in both two cases.

The whole analysis process shows that all the five important parameters we choose will not cause obvious changes in fatty acids but they all have positive effects on the decrease of P. acnes in the whole system.

And surprisingly, for each degree of change, we can see that p3 and K1 both have quite apparent influences on the decrease of fatty acids and P. acnes, which means the system (including engineered bacteria, P. acnes, fatty acid and bacteriocins) represented by the growth simulator model is quite sensitive to these two parameters.


Overall, though we can only display the trend and relative value in the results of our model, the sensitivity analysis is really informative.

Thus, we can improve the effect of our ultimate product by reduce p3 (improving the stability of bacteriocins in the outside environment) and increase k1 (improving the ability of our engineered bacteria to decompose fatty acids). Especially if we can increase k1, the engineered bacteria will be more competitive for nutrient fatty acids, which is reasonable and is exactly what we want to do in the next model.

However, we have to notice that the positive effects resulting from the change of parameters are limited. According to the vertical axis from above figures, the degree of change from 5% to 10% will lead to almost doubling values, but this phenomenon cannot be seen in the change from 10% to 20%. Therefore, as long as we get a satisfying result, there is no need for us to chase the best effect.

3.The flux accelerator model

Part One


Based on the analysis of the The growth simulator model, we propose the design idea of inhibiting the growth of P. acnes by accelerating the decomposition of fatty acids by engineered bacteria. Because of this, we designed the second model - The flux accelerator model. Targeted at improving engineered bacteria’s ability of consuming fatty acids, we build a two-part FBA (flux balance analysis) model[3] based on fatty acid β-Oxidation to see which enzymes should be overexpressed. Fatty acid β-Oxidation is the main process of engineered bacteria decomposing fatty acids, which involves the catalysis of a variety of enzymes. Overexpressing these enzymes will accelerate the decomposition rate of fatty acids, so as to improve the competitive effect of engineered bacteria on fatty acids. However, due to the fact that there is a limitation of metabolic pressure and the system should not be too complicated, it’s not a wise choice to overexpress all enzymes. We establish a flux balance analysis model and draw the conclusion that the enzymes most worthy of overexpression are FadD and FadE.


Apart from the basic assumptions of FBA method, we have identified some additional assumptions based on the characteristics of our project.

1. In any environment, engineered bacteria will finally reach a stable equilibrium under the constraint of stoichiometry.
2. The purpose of this part is to consume fatty acids as much as possible, so exogenic FA (fatty acids) is assumed to be excessive.
3. The rate of each reaction is relatively stable in the Fatty Acid β-Oxidation process.
4. The reaction is complete.

Symbols and definitions

When we construct the metabolic network of Fatty Acid β-Oxidation process, the following symbols are used to represent the reaction process and substances:

Symbols Definitions
M1 acyl-CoA
M2 enoyl-CoA
M3 3-hydroxyacyl-CoA
M4 3-ketoacyl-CoA

With these symbols, Fatty Acids β-oxidation can be expressed as:

Fig7. Pathway of Fatty Acids β-oxidation

Thus, the metabolic network can be described in the following metrological matrix form, with each row representing a substance and each column representing a one-step reaction:

Reaction1 Reaction2 Reaction3 Reaction4 Reaction4* Reaction5
M1 1 -1 0 0 0 1
M2 0 1 -1 0 0 0
M3 0 0 1 -1 -1 0
M4 0 0 0 1 0 -1

The metrology matrix is recorded as:

Set the reaction flux as , each component represents the reaction rate of the corresponding reaction.

Model construction and explanation

System Definition:

Reaction2,3,4,5 form a loop. Reaction1 is the cyclic input substance M1, and reaction4 outputs the substance acetyl coA. In Reaction4, each one-unit participates in the reaction and generates one-unit M4 and one-unit acetyl CoA at the same time. We determine the process of generating acetyl CoA in Reaction4 * as the boundary of reaction flow and Reaction1 as another boundary.


Maximize the consumption speed of exogenic FA, namely, find the maximum value of Reaction1 rate under constraints.

Restraint Conditions:

1. In Reaction4, M4 and acetyl CoA have the same coefficient, which means that v4 = v4*

2. In the final state, the material reaction in the cycle reaches equilibrium, so the variation of M1, M2, M3 and M4 is 0, which means equation constraint:

3. The upper and lower bounds of all reaction Fluxes are shown as follow:

component v1 v2 v3 v4 v4* v5
Upper boundary 1 7.5611 0.8610 0.8610 1.1456 0.7809 0.8610
Upper boundary 2 8.6254 0.7651 0.7651 0.4257 0.6852 0.7651
Upper boundary 3 8.2003 0.7467 0.7467 0.6312 0.0888 0.7467
Lower boundary 0 0 0 0 0 0

Remark: we use the data in the database[4] to obtain the relative data of a unified unit to set our upper and lower bound values, which can save the troublesome steps of measuring accurate values and enhance the operability. The three values in the upper bound represent the data after E. coli is exposed to air for 5 min, 10 min and 15 min. The data is gained from a survey of E. coli transcriptome dynamics during the transition from anaerobic to aerobic conditions.


We use the software MATLAB to solve the optimization problem and compare the values of the target flow vector with the upper bound. It is found that the rate of Reaction2 will improve the restricted target when exposed to the air for a short time, but the target reaction rate is obviously limited by Reaction1 when exposed to the air for 15min. As a result, we should at least increase the expression of FadD, the coenzyme of Reaction1. However, according to the results obtained by previous teams, we found that our results were improvable, so we used the existing flow balance analysis tools to consider β-Oxidation pathway more comprehensively and carefully (Part 2).

Part Two


For the fatty acid β-Oxidation, if the process is only roughly considered as one cycle, many metabolites and metabolic processes will inevitably be ignored. Therefore, in order to simulate the process more accurately, we need to take the whole reaction and metabolite network in E. coli into account and establish a more comprehensive FBA model.

We have obtained a genome-scale metabolic network (GEM) from the Biochemical, Genetic, and Genomic (BiGG) database. iJO1366[5] is a genome-scale reconstruction of the entire metabolic network in Escherichia coli str. K-12 substr. MG1655, a well-researched and well-documented strain that is assumed to be metabolically equivalent to E. coli K-12 TG1. And it accounts for 1366 genes, 2251 metabolic reactions, and 1136 unique metabolites.


Because this model is an improvement of the model before, the basic assumptions are consistent. In addition, there are the following assumptions:

1. Only the fatty acid β-Oxidation is considered.
2. The reaction flux of the decomposition of palmitic acid is used as an index to measure the ability to decompose fatty acids.

Model and explanation

The fatty acid β-Oxidation model is constructed by using the BiGG Model iJO1366 of Escherichia coli str. K-12 substr. MG1655, and the following metabolic network diagram has been obtained:[6][7]

Fig. 8 Map of β-Oxidation model

Aiming at the maximum growth rate of Escherichia coli, we uploaded metabolic data and the following initial metabolic network diagram can be constructed:

Fig. 9 Initial metabolic network


1. In this case, the growth rate of Escherichia coli is 0.982h−1.
2. Color of each reaction:

Min: gray; Med, purple; Max, red, as shown in the Fig. below.

That is, the color represents the relative reaction flux, and the absolute reaction flux is the number on the right side of this step.

3. The direction of the arrow is the defined direction of the forward reaction.
4. Some reaction rates seem to be zero, but in fact, they are due to the number of significant figures, which are actually not zero.


The optimization goal is to maximize the decomposition rate of fatty acids. The 4 results are as follows:

1. Under basic conditions, the maximum consumption of fatty acids is 88.324 mmol/ (gDW*hr). (The amount of the metabolite in millimoles in a gram of bacterial cells (dry weight) over the time of an hour.) That is, it is feasible to increase the consumption of fatty acids under limited conditions.

Fig. 10 The metabolic network diagram when the optimization goal is to maximize the decomposition rate of fatty acids under basic conditions

2. In addition, we set the reaction flux of the specific step as 1 to simulate the situation when FadE is overexpressed. With the maximum decomposition rate of fatty acids as the optimization objective, the results are as follows:

Fig. 11 The metabolic network diagram when the optimization goal is to maximize the decomposition rate of fatty acids and FadE is overexpressed

According to the result, the reaction rate of the decomposition of fatty acid under the circumstances is 86.610 mmol/ (gDW*hr).

3. When we set the reaction fluxes of the reaction that FadB and FadJ are involved in to simulate the overexpression of these two enzymes (set the reaction rate of these two steps as 1), taking the maximum decomposition rate of fatty acid as the optimization objective as well, there is no solution, which indicates that the reaction flux of these two steps cannot be larger if the decomposition rate of fatty acid aims to be the fastest.
The maximum reaction fluxes of these steps are both 1.95*e-14 mmol/ (gDW*hr) if the maximum reaction fluxes of these steps are taken as the optimization objective.

4.In terms of FadA and FadL (setting the reaction flux of these two steps as 1), taking the maximum decomposition rate of fatty acid as the optimization objective, the results are consistent with the FadB and FadJ in result3.

Conclusions of two parts

Based on the above results, the following conclusions can be drawn:

a)From a computational point of view, it is very feasible to increase the consumption of fatty acids.

b)Over-expression of FadD and FadE can improve the ability of our engineered bacteria to decompose fatty acids, while the over-expression in the other two steps of fatty acid β-Oxidation is of little significance and may increase metabolic pressure.

According to the results, only the effect of the most promising enzymes FadD and FadE are overexpressed in our engineered bacteria and tested experimentally.


Due to the limit of time and experimental conditions, our models are not complete to some extent. If time permits, we plan to polish our models and investigate deeper into them in these aspects:

1. Adjust parameters of our growth simulator or add more details to it in accordance with experiment data and try to make it appliable to more conditions that haven’t been considered yet.
2. Learn more about the FBA method and bring up more reasonable and detailed constraints in order to get a better solution.


[2]Momeni, B. , X. Li , and W. Shou . "Lotka-Volterra pairwise modeling fails to capture diverse pairwise microbial interactions." eLife,6,(2017-03-18) 6(2017).
[3]Kenneth J Kauffman; Purusharth Prakash; Jeremy S Edwards (2003). Advances in flux balance analysis. , 14(5), 491–496. doi:10.1016/j.copbio.2003.08.001
[4]Partridge JD, Scott C, Tang Y, Poole RK, Green J. Escherichia coli transcriptome dynamics during the transition from anaerobic to aerobic conditions. J Biol Chem. 2006 Sep 22;281(38):27806-15. doi: 10.1074/jbc.M603450200. Epub 2006 Jul 20. PMID: 16857675.
[5]BiGG Model iJO1366 ( )
[6]Rowe E, Palsson BO, King ZA. Escher-FBA: a web application for interactive flux balance analysis. BMC Syst Biol. 2018 Sep 26; 12(1):84. doi: 10.1186/s12918-018-0607-5. PMID: 30257674; PMCID: PMC6158907.
[7]Heirendt L, Arreckx S, Pfau T, Mendoza SN, Richelle A, Heinken A, et al. Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nature Protocols. 2019 Mar 1; 14(3):639–702.