Dynamic Flux Balance Analysis
Model Description
The previous models characterize the L-cysteine overproduction and toxin–antitoxin kill switch inserts independently. While these inserts are also initially tested separately in the wet lab, the final system is designed to function as a single plasmid within the device. Therefore, developing an integrated model that combines the two modules is essential for understanding the overall behavior of the device over time and within the bioreactor.
The engineered system is inherently dynamic and functions over time. One of the key limitations of the previous models is with constraint-based modeling. The model assumes steady state conditions and the results remain constant, which is not characteristic of the device. Thus, we implemented dFBA to add a time variable to the constraint-based model. Adding a time variable creates the opportunity to link it with the mechanistic/stochastic model to observe how the effects of the two inserts are dependent on one another. A limitation from the mechanistic model was the constant intracellular L-cysteine concentrations used. dFBA addresses the previously mentioned model limitations by linking the different modeling frameworks. The overall aim of the model is to understand how intracellular L-cysteine concentrations from the constraint-based model affect the formation of the transcription factor to activate the toxin and eventually trigger the kill switch from the mechanistic model. This integration enables more accurate prediction of the timing of kill switch activation.
The model was built in Python as an extension of the FBA and all code is deposited in the team’s GitHub repository on the Model Overview page.
Methods
The dFBA extends the results of the FBA by accounting for changes over time and calculating concentrations. Euler’s method was used to implement the dFBA by solving and re-optimizing the FBA model over small-time steps while updating extracellular metabolite concentrations and accounting for nutrient availability.
Euler’s method was implemented through a Python time loop, where the model was optimized using lexicographic optimization and various bounds were updated to set up the subsequent time step.The biomass concentration was calculated using the biomass growth rate predicted by the model at each time step. The overall biomass concentration continually increased at each time step. More specifically, the biomass concentration was modeled to follow different phases of E.coli growth including lag, exponential, stationary, and death phase based on the elapsed time. The biomass growth curve generated by the model is intended to be calibrated based on the growth curves generated experimentally to improve accuracy of the model and L-cysteine export flux predictions.
The values of the concentrations of extracellular metabolites and biomass growth calculated at each time-step were stored in an array to plot and illustrate the behavior of the system over time.
At the end of the time loop, all of the extracellular metabolites including L-cysteine and biomass growth are plotted over time to illustrate the behavior of the system. Furthermore, the dFBA is linked to the mechanistic model through the intracellular L-cysteine accumulation concentrations. By using the same number of time steps and time step duration for the dFBA and the mechanistic model, the intracellular L-cysteine accumulation replaces the previous placeholder constant L-cysteine concentration. The average time to activate the kill switch is recalculated using the updated concentrations to refine the results of the model. We used this integration method to link the separate models together since the independent mechanistic model assumes constant intracellular L-cysteine concentrations. Using the L-cysteine concentrations predicted by dFBA eliminates the need for a separate differential equation, and thus the use of kinetic parameters that are difficult to find and characterize.
Results
The dFBA model without integration with the mechanistic model predicted extracellular metabolites, biomass concentrations, and intracellular/extracellular L-cystiene concentrations over time. However, upon integration with the mechanistic model, the transcription factor concentrations were too low to activate the kill switch. We are currently in the process of troubleshooting the error. The predicted intracellular L-cysteine concentrations from the dFBA reach the target concentration levels of around 300 mg/L, which should be sufficient to activate the kill-switch. The low levels of the transcription factor are likely attributed to the assumption that 8 L-cysteine molecules are required to form the complex and the formulation of the Hill-like approximations used in the mechanistic model. We expect to identify a solution and accurately predict the kill-switch activation time.
Future Directions/Model Calibration
To improve the accuracy of the dFBA, we plan to build and run a calibration program using experimental data. This program will combine a random sampling function and a Mean Squared Error (MSE) cost function. The random sampling function selects values between 10-fold down and 10-fold up of the approximate parameter values in the uncertain parameter space. The uncertain parameter space consists of the modified kcat values and gene expression values in the constraint-based model, the tuning parameters in the CcdR pathway and the expression fluxes in the CcdA/CcdB pathway of the mechanistic model.
The dFBA will then be run using the randomly generated parameters from the random sampling function. The results of the dFBA will then be compared to experimental data, specifically OD600 readings, using an MSE cost function. The MSE cost function calculates the average of the squared differences between the predicted and experimental values. The closer this value is to 0, the more accurate the model predictions are.
After the MSE cost function is run, the sampled parameter spaces will then be analyzed to reconstrain the parameter space. The range that can be randomly sampled will be re-constrained by taking the ranges of parameters of the ten percent parameter spaces with the lowest MSE cost. This process will run iteratively to produce a final optimized parameter space that will be able to accurately represent the results of the benchwork.
Once an accurate model is produced, modularity can be considered. By using literature to find approximate changes gene expression, kcat values, the tuning parameters in the CcdR pathway, and the expression fluxes in the CcdA/CcdB pathway, the behavior of a modified Cysteinator can be approximated. This can then be overlayed to testing feedstocks by changing the medium constraints on the FBA.
In addition to this calibration, we plan to design and develop models to determine how the Cysteinator will interact with the microbes in DF feedstocks, incorporate copy number, and incorporate the rates L-cysteine oxidation to L-cystine.