Modelling Overview

Implementing mathematical and computational modelling methods to understand the co-culture, nisin and PHB


How do we know for a fact that K. xylinus would grow with E. coli or that nisin is able to both be integrated into bacterial cellulose (BC) and still work as an antimicrobial peptide? These were all questions that birthed our team’s three modelling projects in order for us to fully understand and validate certain conditions, factors or systems that needed to be understood before conducting experiments. Modelling also helped us understand conditions that could not be experimentally determined.

Our modelling work consists of three components which includes a co-culture model to describe and understand the workings of our co-culture between K. xylinus, to produce BC and recombinant E. coli to produce either nisin or polyhydroxybutyrate (PHB). The second is a dynamic protein modelling approach that involved both mathematical modelling, docking and molecular dynamics simulations that helped study the production and behavior of our chosen antimicrobial peptide nisin. Lastly, a robust model was modified to help us understand the production and secretion of PHB.

Co-culture Modelling

Co-cultures are complex systems as the interactions between two or more bacteria species causes an increase in the number of factors one has to consider to validate their optimal functioning. This was the reasoning behind creating an exploratory mathematical model that enabled us to look into the dynamics of K. xylinus, and E. coli to allow for the production of BC, PHB and nisin.

The objectives of this model were as follows:

  • Describe the biomass growth of K.xylius and E. coli
  • Understand substrate consumption of the system
  • Describe production of BC and recombinantly produced molecules
  • Understand how other extracellular products affect the growth of both bacteria
  • Compare and contrast monocultures to co-cultures

A System of Ordinary Differential Equations

To start off our modelling , the most important thing to validate was that our bacteria can grow given that they are both dependent on the same source of substrate in this environment. To model the growth, we made use of the growth rate μ as described by Monod's equation (1). Ultimately, this enabled us to derive the following equations:

Equations that describes the biomass growth of K. xylinus and E. coli

The yield coefficient is a constant that can be used to relate cell growth to the rate of substrate consumption. Thus, the yield coefficients of both K. xylinus and E. coli was used to generate an equation to describe the utilization of their substrate, glucose, in the co-culture:

Equation that describes glucose consumption during microbial growth

The next step in modelling was to generate an equation that can help describe the rate and amounts of production of our substances of interest. To start off, we began by modelling BC production. The model assumes that the amount of BC is directly proportional to the amount of K. xylinus present and a product yield coefficient, Y[BC]/X, K. xylinus can be defined (2). Hence, the following equation was derived:

However, there existed difficulties in measuring the amount of K. xylinus in the lab so the wet lab began to assume that BC mass was reflective of how much K. xylinus was present. Since, in some K. xylinus strains, almost 40% of the sole carbon source can be used for BC production, we then modified the model to make the assumption that BC production is dependent on the substrate concentration. So, Monod’s equation was used to implement this modification:

Our protein production model was taken directly from Calleja et al. 2016 which divided the process into three blocks: describing biomass growth, IPTG uptake, and Protein production. The full description on this model can be found on our nisin modelling page.

Lastly, we understood that there might exist other substances secreted by either K. xylinus and E. coli that can affect the environment and the functioning of the other species. We then discovered that K. xylinus produces two acid substances, gluconic and acetic acid, that have the potential of altering the pH of the system which can be detrimental to the growth of the bacteria and limit the production of our desired products. We used the product yield coefficient as a means of relating the production of these acids to K. xylinus biomass. Hence, we derived the following system of equations:

We also derived an equation to relate this to how pH changes over time as these substances are produced.

Simulation Results

The full system of ODEs in the model can be found in the Appendix

In order to solve the system of equations numerically we implemented these equations in Python and made use of the odeint and solve_ivp method from the SciPy package. The python files for all implementations of the model can be found here . The plotted solutions gave rise to the following graphs:

Figure 1. Change in Substrate Consumption for the Co-culture during the simulation run

Figure 2. Change in Bacterial Cellulose in a Co-culture (note that bacterial cellulose amount is used to quantify the amount of K. xylinus)

Figure 3. Change in E. coli biomass in a co-culture

Figure 4. Change in Acetic and Gluconic acid production by K. xylinus in a co-culture

Figure 5. Change in pH of the Co-culture

Explore more on this Model in the Co-culture Modelling page.

Nisin Modelling

nisin production by E. coli has not been thoroughly studied, especially E. coli that are in co-cultures. It is for this reason that we modelled the production of nisin, in order to have an understanding of the production of nisin and how that would be affected when put under the constraints of a co-culture.

Computational tools of docking and molecular dynamics simulations were leveraged to enable us to understand the behaviours and interactions nisin, PHB and BC have with one another. These simulations also helped us study and characterize how nisin’s conformation might change when put under several temperature, pressure and pH conditions.

We implemented a mathematical model to capture nisin production and these were obtained from Calleja et. al 2016 that modeled protein production in fed-batch E. coli cultures. This was an ideal model to use as it implemented our chassis, E. coli, and also considered the fed-batch strategy which is synonymous to the intermittent feeding strategy used in wet lab experiments.

The Calleja model consisted of three components which were the biomass growth mode, the IPTG uptake model and the protein production model (3).

Figure 6. Blocks and Variables of the Model by Calleja et. al (1)

nisin Production

The model takes into account two timestamps; before induction and after induction. During the non-induced stage, only the biomass model is taken into account while after induction, the other models come into play with modifications on the biomass model to account for induction. However, we made modifications to the model to be more reflective of our experimental conditions. One of such is the change in volume equation. The model takes into account a different feeding flow rate equation which the volume is dependent on. As a result, we changed that to reflect our periodic feeding. This was then reflected in all equations that utilize the feeding flow rate. Hence, we present the following equations that describe our nisin production for a monoculture of E. coli:

Non-induced Period:

Induced Period:

We then used the solve_ivp function in Python to solve the system of equations in order to get plots that can help visualize the changes in the modelled parameters.

Figure 7. Change in E. coli biomass during the simulation run

Figure 8. Change in Substrate Consumption during the simulation run

Figure 9. Protein production during the simulation run

Non-induced period of the simulation is shown in blue and induced periods of the simulation is shown in red.

Although initially produced in a monoculture, E. coli would be in a co-culture with K. xylinus. Therefore, it was important that we understood how this could affect the production of nisin. So we added changes into the model for a co-culture and plotted the model results. This resulted in the following:

Figure 10. Change in Substrate Consumption forthe Co-culture during the simulation run

Figure 11. Change in Bacterial Cellulose in a Co-culture (note that bacterial cellulose amount is used to quantify the amount of K. xylinus)

Figure 12. Change in E. coli biomass in a co-culture

Figure 13. Change in Acetic acid production by K. xylinus in a co-culture

Figure 14. Change in Gluconic acid production by K. xylinus in a co-culture

Figure 15. Protein production in the co-culture during the simulation run

Non-induced period of the simulation is shown in blue and induced periods of the simulation is shown in red.

Explore more on this Model in the nisin Modelling page.


nisin’s placement in Cellucoat is solely dependent on the interactions nisin has with BC. To look into these interactions and confirm that they are sufficient to immobilize nisin unto BC, we performed docking simulations. However, PHB, the bioplastic added to our BC to give it additional mechanical properties, would also come in contact with nisin in the same environment. This was another interaction we studied through docking simulations.

Our docking simulations were conducted on a web server CB-DOCK 2 (4) that is used for protein-ligand docking. Given the three-dimensional structure of a protein and ligand, it integrates cavity detection, a docking algorithm and homologous template fitting to predict the sites and affinity that a ligand has to a protein (4).

The best out of five possible binding conformations was taken and used as the basis of our analysis. The Vina scores were then used as our metric of measuring binding affinity. The score works in such a way that the more negative it is, the less energy that is needed to initiate binding and the better the binding affinity.

Protein + Ligand Vina Score Cavity Volume Center Docking Size
nisin + BC -5.4 83 13, 16, -12 33, 33, 33
nisin + PHB -4.1 83 13, 16, -12 27, 27, 27

Figure 16. 3D Visualization of nisin binding unto a chain of BC (Contact residues of nisin to BC shown in red and other parts in green. BC chain shown in red and white)

Figure 17. 3D Visualization of nisin binding unto a chain of PHB (Contact residues of nisin to PHB shown in red and other parts in green. PHB chains shown in yellow and red)

After conducting docking, we analysed the contact residues of PHB and BC with nisin and discovered that the best binding conformation for both is at the same amino acid residues. However, BC (-5.4) has a lower Vina score than PHB (-4.1) meaning that it has a higher binding affinity. So empirically, nisin would favourably bind more to BC than PHB. Thus, allowing for immobilization of nisin into BC.

In nisin’s antimicrobial mechanism, when nisin reaches the cell membrane surface, rings A and B at the N-terminus of nisin bind to the pyrophosphate group of lipid II. The C-terminal region inserts into the membrane to form a pore using a flexible mechanism. The aforementioned regions and rings do not overlap with the binding sites of PHB and BC.

Molecular Dynamics

In the production-to-use pipeline of Cellucoat, it is going to be placed under certain environmental conditions that derail from regular and optimal conditions. nisin’s functionality under these conditions could not be validated experimentally due to time and resource constraints. So we turned to conducting molecular dynamics simulations to validate this.

Two such conditions was that nisin was going to be put under the high temperature and pressure conditions of the autoclave or other sterilizing machines after the production of BC so as to kill off any bacteria present and lyse open E. coli to release nisin as found on the Co-culture page. On the other end of the spectrum, we understand that fruits when in use and in storage are often placed in cold conditions in fridges and refrigerators.


To understand how this variance in conditions would affect nisin, we conducted three molecular dynamic simulations. We assumed in our simulation that the only parameters that change from the control when under other conditions are temperature and pressure. The first served as a control that modelled nisin in room temperature and under normal atmospheric pressure conditions. The second was simulating nisin when under autoclave temperature and pressure conditions and the last simulation involved nisin in fridge conditions.

To account for these temperature and pressure changes, .mdp files in the simulation run which contain parameters that control the pressure and temperature of the simulation were changed. They were changed as follows:

Simulation Simulation Temperature Simulation Pressure
Normal Conditions (control) 27 ° C or 300K 14.5 psi or 1 bar
Autoclave Conditions 0 ° or 273K 15 psi or 1.03421 bar
Fridge Conditions 0 ° or 273K 28 psi or 1.93053 bar

Table 2. Temperature and Pressure values of the several molecular dynamics simulations on nisin (5 - 8).


Figure 18. 3D Visualization of the molecular dynamics simulation of nisin under normal conditions

Figure 19. 3D Visualization of the molecular dynamics simulation of nisin under autoclave conditions

Figure 20. 3D Visualization of the molecular dynamics simulation of nisin under fridge conditions

After conducting the simulations we wanted to analyse how the simulation condition affects nisin’s conformation. So to do this we expressed the change in conformation terms of how the distance between key amino acids are changing over the period of the simulation. Two amino acids we studied were cysteine 7 and serine 33 which are the first and last amino acid residues in nisin that are in contact with BC as resulted from the docking simulation conducted.

So we used a software called VMD to generate the distance between the two amino acids over the frames of the simulation. These were extracted as a file of the distance over frames and also a histogram of the distance.

The results are as follows:

Figure 21. Distance between Cysteine 7 and Serine 33 residues in nisin under normal conditions

Figure 22. Distance between Cysteine 7 and Serine 33 residues in nisin under fridge conditions

Figure 23. Distance between Cysteine 7 and Serine 33 residues in nisin under fridge conditions

Figure 24. Distance distribution of nisin under normal, fridge and autoclave conditions between Cysteine 7 and Serine 33 residues

Again it was also important to validate that nisin is still able to bind unto BC so that it can be adequately integrated in Cellucoat. As we had done before, we conducted docking simulations with the nisin structure that had undergone the molecular dynamics simulation. We removed all the water molecules and ions and got the protein structure of nisin and conducted docking with it to BC. These resulted in the following:

Simulation Condition Vina Score Cavity Volume Center Docking Size
Autoclave (nisin + BC) -5.3 76 40, 40, 8 33, 33, 33
Fridge (nisin + BC) -4.9 42 30, 25, 15 33, 33, 33

Table 3. Docking Results of nisin and BC under autoclave and fridge conditions

From the docking simulations, the reductions in binding efficacy of nisin and BC were below 5%. Thus, we can confirm through the model that nisin would still be able to bind unto BC and still have similar conformations to itself when under normal conditions while in the conditions of the autoclave or fridge.

Explore more on this Model in the nisin Modelling page.

PHB Modelling

PHB is produced in E. coli to be incorporated into BC to allow for it to attain additional mechanical properties that make it a more comparable alternative to traditional plastic packaging. This is because on its own, BC is not mechanically sufficient to be used in the produce packaging industry.

Our production of PHB leverage iGEM Calgary 2017’s part (BBa_K934001) which created a PHB production and secretion system. It included a process whereby Hly-A tagged phasin is produced in E. coli, allowing for the secretion of PHB via the type I secretion pathway (9). Phasin has a role in allowing for PHB granule size regulation and the ability of the cell to secrete PHB granules. So, our work involved proposing new part designs that leveraged the changing of the strength of the ribosomal binding site (RBS) of the phasin coding sequence. This is done to have an effect on phasin expression levels, which would allow for more phasin production and increase the secretion of PHB into the BC. Ultimately, this would imply a change in the mechanical properties of the material, making it a more suitable alternative in the packaging industry.

To have an understanding of PHB production and the phasin secretion system, although this is possible by conducting lab experiments, we wanted to use modelling to validate this system.

Modelling Method

Dixon created a predictive model to describe and ultimately optimize the metabolic pathways involved in the production of PHB in E. coli (10). The model is implemented in MATLAB 2010 using the simbiology model builder module. The model consists of the main metabolic pathways used in PHB production which are glycolysis, acetyl coenzyme A (acetyl-CoA) synthesis, tricarboxylic acid (TCA) cycle, glyoxylate bypass, and PHB synthesis (10). The reactions in these pathways are mathematically described using kinetic mechanisms of enzymes that are in charge of catalysing the reactions. These reactions are dependent on enzyme concentration which is solely a result of the rate at which the enzymes are transcribed, translated and degraded. Thus, the model also took into account how these enzymes are being produced as a function of the transcription factors and promoters in charge of the enzyme’s production.

Modifications to Reflect our Design

We modified the model to capture all the details we needed in adequately describing our PHB production and phasin secretion system. One of such changes we made was adding in a component that accounted for phasin production in the cell. The original model accounts for enzyme production as a function of promoters and transcription factors. PHB is produced as a result of three enzymes which are β-ketothiolase, acetoacetyl-CoA reductase and poly-β-hydroxybutyrate polymerase (PHB synthase) that convert acetylCoA to PHB. The genes that encode for these enzymes are transcribed using the lac promoter and the model states that the transcription factors Crp, H-NS, and LacI are used in regulating the lac promoter. Furthermore, this translates to the rate and amount of PHB synthesis.

Since phasin is inserted on the same plasmid as these enzymes in charge of PHB production, we assumed that the same transcription factors would also be used in regulating its transcription. Thus, we based phasin production on these transcription factors as is the manner of the model.

Equation 1. Equation 1 describes the concentration of protein as a function of multiple transscription factors (2). This was used to describe the production of phasin and the other enzymes involved in PHB synthesis.

After incorporating phasin production into the model, it was necessary that we translate its production to how it interacts with PHB in order to regulate its size and allow for its secretion from E. coli. We assumed that PHB and phasin formed a complex which its production could then be related to a kinetic mechanism. Using mass action law, we then created an equation that we used to describe the reaction of PHB and phasin to form the PHB-phasin complex.

Reaction 1. Reaction 1 describes the binding of intracellular PHB and Phasin to form a "PHB-Phasin complex" that can be secreted using type 1 secretion system in E. coli

As expanded upon earlier, the incorporation of phasin allows for PHB granules to be secreted out of the cell and we wanted to capture this in the model. The original model assumes PHB stays inside the E. coli cell, not secreted. So we made a new component and added it to model a reaction relating to secretion of the PHB-phasin complex and the breaking down of the complex into PHB granules and phasin.

Reaction 2. Reaction 2 describes "PHB-Phasin complex" being seperated into extracellular PHB and Phasin(Pout)

Lastly, an additional change we made to the model was to account for the degradation of both PHB and phasin, as the model had not done so earlier although it had captured that of the enzymes. This was to ensure that when simulations were run, only the net amount was revealed, hence making the model a more accurate depiction of the actual system.

Simulation results

After implementing all our changes and changing parameters for some of the other reactions, we ran the resulting model as a simulation. The resulting plots are as follows:

Figure 25. Production of β-ketothiolase (PhaA), Acetoacetyl-CoA (PhaB) reductase and Poly-β-hydroxybutyrate polymerase (PHB synthase - (PhaC)) that convert AcetylCoA to PHB

Figure 26. Change in PHB concentration

Figure 27. Change in intracellular Phasin (Pin) concentration

Figure 28. Change in "PHB-Phasin" complex

Figure 29. Change in extracellular PHB (PHBout) concentration

Analysing the simulation run gave us an understanding that the Phasin secretion system is a sufficient system that can lead to the secretion of PHB. One notable key thing is that at the 150,000 seconds mark of the model it is evident that the intracellular concentration of PHB is reduced to zero. This resulted in the decrease in the production “PHB-Phasin” complex. This can lead to a reduction in the amount of PHB secreted and can decrease the desired mechanical properties of BC. Hence, there has to be a mitigation strategy implemented to ensure that the production of PHB is not halted. This can be such as ensuring increased amounts of substrate being given to the chassis during fermentation. In the wet lab the implementation of an intermittent feeding strategy achieves just this. This is because K. xylinus can grow for long periods and to sustain a consistent production of PHB, there has to be sufficient amounts of substrate.

Explore more on this Model in the PHB Modelling page.

Find all equations and parameters used in the models here


  1. Monod J. The growth of bacterial cultures. Annu Rev Microbiol [Internet]. 1949 Oct [cited 2022 Sep 19];3(1):371–94. Available from:
  2. Notebook [Internet]. [cited 2022 Sep 19]. Available from:
  3. Calleja D, Kavanagh J, de Mas C, López-Santín J. Simulation and prediction of protein production in fed-batch E. coli cultures: An engineering approach: Simulation of E. coli Protein Production. Biotechnol Bioeng [Internet]. 2016 Apr [cited 2022 Oct 10];113(4):772–82. Available from:
  4. Liu Y, Yang X, Gan J, Chen S, Xiao ZX, Cao Y. CB-Dock2: improved protein-ligand blind docking by integrating cavity detection, docking and homologous template fitting. Nucleic Acids Res. 2022 May 24;gkac394.
  5. Autoclave overview [Internet]. [cited 2022 Oct 10]. Available from:
  6. Fruits and vegetables - optimal storage conditions [Internet]. [cited 2022 Oct 10]. Available from:
  7. Freon™ refrigerants for residential refrigeration [Internet]. [cited 2022 Oct 10]. Available from:
  8. Refrigerant pressures, states, and conditions | achr news [Internet]. [cited 2022 Oct 10]. Available from:
  9. Part:BBa K2260002:Experience - [Internet]. [cited 2022 Sep 11]. Available from:
  10. Dixon A. Designing predictive mathematical models for the metabolic pathways associated with polyhydroxybutyrate synthesis in escherichia coli. 2011 [cited 2022 Oct 10]; Available from: