MATLAB code for modelling how synthetic gene circuits interact with the host cell's (E. coli) native genes and affect its growth rate, used in the manuscript 'A coarse-grained bacterial cell model for resource-aware analysis and design of synthetic gene circuits'.
The folders are organised as follows:
Contains scripts implementing the model of an E.coli cell. These include:
- cell_simulator.m - a Matlab class object enabling simulations of the host cell. The expression of synthetic circuits can be simulated by loading the 'heterologous gene' and 'external input' modules (see het_modules and ext_inputs folders). Note that the associated modules' parameters must be pushed into the main simulator's memory using the function push_het() every time they are altered.
- cell_params.m - provides default values of all parameters describing the host cell
- cell_init_conds.m - provides the default initial conditions for simulating the host cell model
- cell_formulae.m - contains a collection of formulae for rate and activation functions used by the cell simulator
- get_steady.m - an auxiliary function that allows to determine the cell's steady state in given conditions by running cell_simulator simulations.
Contains scripts defining class objects that allow to model the expression of different heterologous gene circuits. These include:
- no_het.m - no heterologous gene being expressed. Due to being 'empty', this script can be copied and filled in to describe a synthetic gene circuit of interest
- one_constit.m - one constitutive heterologous gene
- two_constit.m - two constitutive heterologous genes
- t7_selfact.m - a non-cooperative self-activator circuit consisting of a T7 RNA polymerase that transcribes its own gene
- two_switches.m - two self-activating heterologous genes acting as bistable switches. If the genes impose a sufficiently high burden on the cell, winner-takes-all behaviour means that the activation of one switch prevent the activation of the other
- pi_controller.m - a proportional-integral feedback controller for maintaining a constant extent of competition for ribosomes in the cell
Contains scripts defining class objects that describe external inputs administered to the cells, such as chemical inducer concentrations or light stimuli. These include:
- no_ext.m - no external signal. Due to being 'empty', this script can be copied and filled in to describe an externa input of interest
- constant_inducer.m - constant concentration of a chemical inducer
- step_inducer.m - a step increase in the concentration of a chemical inducer
- pulse_inducer.m - time-limited step pulse above or below the baseline concentration of a chemical inducer
Running the corresponding scripts allows to reproduce the simulation results presented in the publication. Files not named fig(...).m are auxiliary scripts required by some of the figure-generating scripts.
Contains scripts allowing to fit the model's parameters to experimental data obtained by Scott et al. [^1] and processed as descroibed by Chure et al. [^2] using the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm [^3]. Note that these scripts require the DREAM package to work.
- dream_fit.m - run the DREAM algorithm and record the outcome in a .mat file
- dream_model.m - defines the function for the model's log-likelihood given the experimental data for an input set of parameter values. This function is called while sampling the chains during DREAM fitting
- dream_modelfun.m - defines the function for finding the sum of squared errors between the experimental measurements and model predictions for an input set of parameter values. This is required to calculate the log-likelihood
- DREAM.mat - stores the sampled chains obtained as outcome of running the the algorithm for 20,000 steps in total (10 parallel chains tun for 1,000 steps each)
- DREAM_postproc_amended.mat - script for postprocessing the fitting outcome (find the mode of the sampled distribution, plot the sampled chains, etc.). This script has been taken from the original DREAM package, but minor amendments have been made to it so as to improve the output plots' interpretability
- update_a_a.mat - script for updating the value of the a_a parameter (assigned an order-of-magnitude estimated value for the duration of MCMC fitting) based on the fitting outcome
Experimental data used to fit the parameters and compare model predictions with real-life measurements, taken from [^1] and [^2]. The text file annotation.txt explains the meaning of each dataset present.
The code was run and ensured to be functional with Matlab R2022a and Matlab R2023b on PCs running on Windows 10 Pro 22H2, Windows 11 Home 22H2 and Ubuntu 20.04.6.
No non-standard hardware is required to run the code. Additional software requirements include the Matlab Parallel Computing Toolbox 7.6 and the DREAM Matlab package, part of the HydroSight toolbox (https://github.com/peterson-tim-j/HydroSight; date last accessed: 10 September 2022).
To install and run the code, download this repository as a zip file, extract its contents into a folder, and open it using Matlab. Provided that all other required software has been installed, the process should take less than 5 minutes on a typical PC.
To reproduce the results displayed in the figures from the article and its supplementary information, run the script with the corresponding name in the 'figures' folder (e.g. fig3a.m will reproduce the main text's Figure 3a, figS1.m will reproduce Supplementary Figure 1...). The obtained plots should match the figures in the text.
Different scripts will have different runtimes on a typical PC. However, with the exception of the scripts listed below, they should not take more than 10 minutes.
-
dream_fit.m - parameter fitting; takes
$\approx 48$ hours provided that the computer can run 10 parallel pool workers at once. Raw parameter fitting outcome data given in DREAM.mat -
fig9.m - demonstration of how the integral controller proposed in the article enforces modularity of synthetic biology designs; takes
$\approx 30$ minutes. -
figS2a.m, figS7abc.m - scripts that invove sampling parameters from the fitted distribution and simulating the model with them; take
$\approx 20$ minutes. -
figS2e.m - model parameter sensitivity analysis; takes
$\approx 3$ hours provided that the computer can run 12 parallel pool workers at once. -
figS5defg.m - bifurcation analysis of a non-cooperative self-activator; takes
$\approx 2$ hours. -
figS7def.m - one trajectory takes
$72-96$ hours to simulate. Multiple trajectories can be sampled at the same time using the Matlab Parallel Computing Toolboox.
For simulating the expression of an arbitrary heterologous gene circuit by the host cell, define the circuit model by filling the template file no_het.m in the 'heterologous modules' folder according to the instructions in the template script's comment lines. If any external input (e.g. inducer pulse) is provided to the system, it can be defined by filling the template file no_ext.m in the 'external inputs' folder according to the instructions in the template script's comment lines -- otherwise, the unmodified no_ext.m script should be used. The system's behaviour can then be simulated by running:
sim=cell_simulator; % initialise simulator
sim=sim.load_heterologous_and_external('heterologous_module_name','external_input_module_name'); % load heterologous gene expression module
sim.tf = simulation_time; % simulation will be run from 0 to simulation_time hours
sim=sim.simulate_model(); % simulate
REFERENCES
[^1] Matthew Scott, Carl W. Gunderson, Eduard M. Mateescu, Zhongge Zhang, and Terence Hwa. Interdependence of cell growth and gene expression: Origins and consequences. Science, 330(6007):1099–1102, 2010.
[^2] Griffin Chure and Jonas Cremer. An optimal regulation of fluxes dictates microbial growth in and out733 of steady-state. eLife 12, e84878, 2023.
[^3] Jasper Vrugt. Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation. Environmental Modelling & Software. 75:273-316, 2016.