Reasoning

When deciding how to model our kill switch, we were initially unsure of how to proceed. We began by defining our goals for what a useful model of our system would be able to provide.

• How long on average does the killswitch take to activate?
• What is the average internal concentration of cysteine at that time?
• What is the desired promoter of ccdA such that we ultimately reach our target concentration at cell death?

These 3 questions were our guide in ultimately determining how our model should function and what method we would use to ultimately understand this system. I would first like to note that the first 2 questions are questions of averages, implying that population dynamics for the model are a must. As a result, we also knew it was necessary that our model had to have some form of induced variability; otherwise, we would achieve the same basic result each time and receive limited information. For this reason, we decided to design a hybrid model, one that uses both a deterministic and stochastic approach to ultimately characterize our system. We knew our best option for the deterministic portion would be Michaelis-Menten kinetics, though we were momentarily stumped on how to introduce variability into the system to increase the accuracy of our model.

Assumptions

Our examinations of the individual parts of the system led us to recognize a very important fact about our transcription factor for ccdB, that being that there isn't a lot of it in the system. This leads us to our first assumption.

• The expression of ccdB is stochastic

The justification for this stems from our current understanding of the CcdR, CcdA* (the cyteine induced on CcdB) system, where our transcription factor is an octamer consisting of 8 ccdR molecules and 2 cysteine molecules (Gao et. al 2021)[1]. Couple this with a low constitutive expression of ccdR, and there are likely on the order \(10^1\) molecules in any given cell. This fact implies that we cannot assume that at any given time, ccdB is being expressed, and therefore, we must model this process stochastically. However, it is unfair to justify this by intuition alone, and we have found that our original assumption was correct, and this can be seen in Section 4. For our next assumption, we are again looking at the ccdR system, and we find that this system is not well characterized. The little information we have on the system comes from Gao et. al, but it does not provide any solid values for the many parts of the system, including, most importantly, formation and degradation fluxes. For this reason, we knew Michaelis-Menten kinetics would not be applicable; however, there is a commonly used technique for scenarios of this type, which involves using a modified Hill equation to characterize the system. So assumption 2 is

• The ccdR system can reasonably be characterized by modified Hill equations.

Much of our justification for this approach stems from existing usages in literature. While many papers have employed this technique, a couple of which we have found useful are "A General Model for Toxin-Antitoxin Module Dynamics Can Explain Persister Cell Formation in E. coli" (Galens et al., 2012)[2] and "A synthetic oscillatory network of transcriptional regulators (Elowitz and Liebler, 2000)[3]. We do, however, recognize the limitations of this method in that it can often overestimate our desired values. For this reason, we have added unitless tuning parameters that allow us to tune our model more effectively to our experimental results. These equations will be discussed later when we talk about our methods. The third and final major guiding assumption for our mechanistic model is

• Cysteine production and therefore accumulation is uniform throughout the population

Without this assumption, we would have to introduce variability into our dFBA to account for individual cells expressing at different rates based on localized macromolecule concentrations. In this assumption, we also assume that our environment is well-mixed and that uptake rates across the population are generally the same. This enables us to apply the same cysteine production and accumulation data to the entire population.

Methods

Below, we discuss in detail the methods by which we have constructed this model, broken up into each section for convenience. When a minor assumption is made, a small justification will be provided unless already provided in the assumption section. Each section will begin with the relevant equations, which will then be discussed in detail, both in their derivation and their relevance to the broader model.

CcdR Pathway

Relevant Equations: \[ \frac{dR_{mRNA}}{dt}=k_{RmRNA}-k_{degM}R_{mRNA} \] \[ \frac{dR}{dt}=k_{Rtrans}R_{mRNA}-k_{Rdeg}R-\frac{2R^3}{K_{d2}+R^2} \] \[ \frac{dR_{2}}{dt}=k_{2form}\frac{R^3}{K_{d2}+R^2}-k_{2use}\frac{2R_{2}^3}{K_{d4}+R_{2}^2} \] \[ \frac{dR_{4}}{dt}=k_{4form}\frac{R_2^3}{K_{d4}+R_2^2}-k_{4use}\frac{2R_{4}^3}{K_{d8}+R_{4}^2} \] \[ \frac{dR_{8}}{dt} = k_{8form}\frac{R_{4}^3}{K_{d8}+R_{4}^2}-k_{8use}\frac{(R_{8}[Cys]^4)^2}{K_{dTF}+(R_{8}[Cys]^8)} \] \[ \frac{dR_{TF}}{dt}= k_{TFform}\frac{(R_{8}[Cys]^4)^2}{K_{dTF}+(R_{8}[Cys]^8)} \] Necessary Assumptions:

• Oligomerization is treated using Hill-like approximations to capture saturation/cooperativity without detailed kinetics.
• \(k_{use}\) and \(k_{form}\) parameters scale the effective rate; they don't represent actual kinetic constants unless calibrated.
• Stoichiometry is explicitly conserved in flux terms.
• All components are assumed well-mixed and cytoplasmic unless specified otherwise.
• Binding of all substrates is perfectly cooperative

Parameter	Value	Definition	Source
\(k_{RmRNA}\)	\(1 \times 10^{-14}\)	Formation rate of CcdR mRNA	Guess based on tuned CcdA parameter
\(k_{degM}\)	0.0203	Degradation rate of CcdR mRNA	From Gelens et al.
\(k_{Rtrans}\)	0.09	Translation rate of CcdR	Guess
\(k_{Rdeg}\)	0.00115524	Degradation rate of CcdR protein	Guess based on Gelens et al.
\(k_{d2}\)	\(8.4 \times 10^{-18}\)	Dimer dissociation constant (\(K_d\))	From Prodigy MD simulation
\(k_{2f}\)	1	Dimerization tuning parameter	Tuning parameter
\(k_{2u}\)	1	Dimer usage tuning parameter	Tuning parameter
\(k_{d4}\)	\(8.9 \times 10^{-9}\)	Tetramer dissociation constant (\(K_d\))	From Prodigy MD simulation
\(k_{4f}\)	1	Tetramerization tuning parameter	Tuning parameter
\(k_{4u}\)	1	Tetramer usage tuning parameter	Tuning parameter
\(k_{d8}\)	\(7.9 \times 10^{-16}\)	Octamer dissociation constant (\(K_d\))	From Prodigy MD simulation
\(k_{8f}\)	1	Octamerization tuning parameter	Tuning parameter
\(k_{8u}\)	1	Octamer usage tuning parameter	Tuning parameter
\(k_{dTF}\)	\(1 \times 10^{-15}\)	Dissociation constant for cysteine binding	From Prodigy MD simulation
\(k_{TFf}\)	1	Cysteine–TF binding tuning parameter	Tuning parameter
\(O\)	cyst	Intracellular cysteine concentration	Input from the dFBA

Chart containing parameters and their sources for the CcdR system (WIP)

Let's start by discussing the derivation of our equations. We first start with the Hill equation. \[ \theta=\frac{[L]^n}{K_d+[L]^n} \] Where L is the concentration of the ligand, \(K_d\) is the dissociation constant associated with the system, and n is the Hill coefficient, which is representative of the level of cooperativity between the protein and ligand. \(\theta\equiv\frac{[P]_{bound}}{[P]_{total}}\), therefore equation 7 can be rewritten as \[ [P]_{bound}=\frac{[L]^n[P]}{K_d+[L]^n} \] Additionally, because we are assuming perfect cooperativity, the Hill coefficient is just equivalent to the number of molecules that are present within the reaction. Therefore, in the case of the dimerization for ccdR, we say n=2 to illustrate perfect binding between molecules. \[ r_{f2}(R)=\frac{R^3}{K_{d2}+R^2} \] We then subtract out the use of the dimer in the formation of the tetramer with a coefficient directly preserving the stoichiometry. \[ r_{f2}(R)-r_{u2}(R_2)=\frac{R^3}{K_{d2}+R^2}-\frac{2R_{2}^3}{K_{d4}+R_{2}^2} \] From here, we then add in both form and use coefficients each with unit \(s^{-1}\), giving u.s \[ \frac{dR_{2}}{dt}=k_{2form}\frac{R^3}{K_{d2}+R^2}-k_{2use}\frac{2R_{2}^3}{K_{d4}+R_{2}^2} \] This same derivation is then repeated for our other coupled ODEs. The above system is the first acting portion of our model, as its outputs are necessary for the formation of our Markov chain for the expression of CcdB. It takes the output of internal cysteine from our dFBA and uses it to calculate our final concentration of transcription factor, which we discuss in our results section.

Stochastic Model

Continuous Probabilistic Equations \[ \partial_tP(c_1,t|x,t_o)=-\lambda_1P(c_1,t|x,t_o)+\lambda_2P(c_2,t|x,t_o) \] \[ \partial_tP(c_2,t|x,t_o)=\lambda_1P(c_1,t|x,t_o)-\lambda_2P(c_2,t|x,t_o) \] \[ \vec{P}=\begin{bmatrix} \partial_tP(c_1,t|x,t_o) & \partial_tP(c_2,t|x,t_o) \end{bmatrix} \] \[ W=\begin{pmatrix} -\lambda_1 & \lambda_2 \\ \lambda_1 & -\lambda_2 \end{pmatrix} \] \[ \frac{d\vec{P}}{dt}=W\vec{P} \] Discretized Probabilistic Equations \[ T_{bind}=1-e^{-k_{on}[TF]\Delta t} \] \[ T_{unbind}=1-e^{-k_{off}\Delta t} \] \[ M=\begin{pmatrix} 1-T_{bind} & T_{bind} \\ T_{unbind} & 1-T_{unbind} \end{pmatrix} \] Necessary Assumptions:

• The probability of binding and unbinding is well modeled by Poissonian statistics
• Our rate constants are of the form \(k_{on}[TF]\) and \(k_{off}\) repsectively
• Our process is Markovian and therefore the transition is only governed by the current state

Parameter	Value	Definition	Source
\(k_{on,TFP}\)	1×108	Binding rate constant for transcription factor–promoter interaction	Estimated typical protein–DNA binding on-rate
\(k_{off,TFP}\)	0.01	Dissociation rate constant for transcription factor–promoter interaction	Estimated by magnitude for a typical DNA-protein dissociation rate

Transcription factor–promoter binding and unbinding rate constants.

Let's begin with the derivation as necessary. We begin by deriving equation 16, though not in vector form \[ \frac{dP}{dt}=-kP(t) \] where k is our rate constant and P(t) is the cumulative probability of some event not happening, we then use the semigroup property of Markovian systems, which form a semigroup under multiplication, such that Letting M(t) be the transition matrix describing the time evolution of the system for continuous t \[ M(t+\Delta t)=M(t)M(\Delta t) \] Therefore \[ P(t+\Delta t)=P(t)P(\Delta t) \] Or in words the probability that the event does not happen up to some \(t+\Delta t\) is equivalent to the probability it doesn't happen up to t times the probability it doesn't happen over some additional discrete time. Additionally, over this discrete time step we know that \[ P(\Delta t)=1-k\Delta t \] where k is the rate of the functional decay. By substituting this expression into Eq. 22 we find that \[ P(t+\Delta t)=P(t)(1-k\Delta t) \] By rearranging terms and taking the limit as \(\Delta t\) approaches 0, we recover Eq. 20. By solving the ODE in Eq. 20 we retrieve the result for P(t) which is as follows \[ P(t)=e^{-kt} \] Which can then be discretized by setting \(t=\Delta t\). To solve for our rate constant,s we perform dimensional analysis. For the binding probability, we know the exponential term must be unitless; therefore, k must have units 1/s. We also know that the probability should increase as a function of the concentration of our transcription factor and as a function of the \(k_{on}\) value for the reaction. \(k_{on}\) has units \(M^{-1}s^{-1}\), therefore our final binding probability must be \[ T_{bind}=1-e^{-k_{on}[TF]\Delta t} \] and by a similar argument, the unbinding probability must be \[ T_{unbind}=1-e^{-k_{off}\Delta t} \] These probabilities can then be placed in the broader transition matrix M which will be used later in our coded 2 state Markov chain where \[ M=\begin{pmatrix} 1-T_{bind} & T_{bind} \\ T_{unbind} & 1-T_{unbind} \end{pmatrix} \] The Markov chain is ultimately what governs the expression of ccdB against ccdA, which under normal expression would generally outweigh ccdA greatly. However, due to this stochastic expression, we see that it is truly governed by the random fluctuation of states within this system. It can be seen that the probability of unbinding is constant, though the probability of binding increases with the concentration of our transcription factor, which is affected by the accumulated amount of intracellular cysteine. This goes to show that the accumulation of intracellular cysteine is what directly governs the activation of our kill-switch and that it is not reliant on extraneous factors.

CcdA and CcdB Pathway

Relevant Equations \[ \frac{dA_{mRNA}}{dt}=k_{AmRNA}-k_{degM}A_{mRNA} \] \[ \frac{dB_{mRNA}}{dt}=k_{BmRNA}-k_{degM}B_{mRNA} \] \[ \frac{dA}{dt} = k_{\text{Atrans}} A_{\text{MRNA}} - k_{\text{Adeg}} A - k_{\text{bind}} A B + 2F k_{\text{Adeg}} C_2 + F k_{\text{Adeg}} C + k_{\text{Cdeg}} C \] \[ \frac{dB}{dt} = k_{\text{Btrans}} B_{\text{MRNA}} - k_{\text{Bdeg}} B - k_{\text{bind}} A B + k_{\text{Cdeg}} C \] \[ \frac{dC}{dt} = k_{\text{bind}} A B + k_{\text{Cdeg}} C_2 - k_{\text{Cdeg}} C - F k_{\text{Adeg}} C - k_{\text{bind}} C B - k_{\text{Bdeg}} C \] \[ \frac{dC_2}{dt} = k_{\text{bind}} B C - k_{\text{Cdeg}} C_2 - F k_{\text{Adeg}} C_2 - k_{\text{Bdeg}} C_2 \] Necessary Assumptions:

• We assume that Michaelis-Menten Kinetics is an effective model for our system
• We assume CcdA expression has reached steady state, and CcdB expression is minimal as compared to when the cysteine over-expression pathway is active
• As with any system consisting of Michaelis-Menten kinetics we assume the system is acting under the Quasi-steady state regime

Parameter	Value	Definition	Source
\(F\)	0.2	Percentage of decay rate for antitoxin within the complex	Gelens et al.
\(k_{Cdeg}\)	0.00000714972	Rate of complex degradation	Gelens et al.
\(k_{bind}\)	0.055	Rate of complex formation	Gelens et al.
\(k_{AMRNA}\)	0.033012	Formation rate of CcdA mRNA	Guesstimation; real value ≈ 0.02751–0.044016
\(k_{degM}\)	0.00203	mRNA degradation rate	Gelens et al.
\(k_{Adeg}\)	0.00115524	Degradation rate of CcdA	Gelens et al.
\(k_{Atrans}\)	0.139	Translation rate of CcdA mRNA	Gelens et al.
\(k_{BMRNA}\)	variable (\(v_{kBMRNA}\))	CcdB mRNA formation rate (depends on transcriptional state)	Guesstimation
\(k_{Btrans}\)	0.033	CcdB translation rate	Gelens et al.
\(k_{Bdeg}\)	0.00577623	Degradation rate of CcdB	Uniprot

Chart containing parameters and their sources for the CcdA/CcdB system.

Michaelis-Menten kinetics is an extremely common approach that draws heavily from experimental data for rate constants, which is why it was the perfect fit for characterizing the CcdA/B system, as it is already well characterized. As such, much of this information was drawn from the 2020 iGEM team from BITS Pilani. Below is a short overview of Michaelis-Menten kinetics and how we applied it to our system.

Formulation of the Michaelis--Menten Framework in the CcdA/B System

The system of equations describing the CcdA/CcdB toxin--antitoxin module is constructed under the assumption that the underlying biochemical interactions can be effectively modeled using Michaelis--Menten kinetics. In this framework, reaction rates are governed by enzyme--substrate--like dynamics, where complex formation and dissociation occur much faster than the synthesis and degradation of the species involved. This allows the use of quasi--steady state approximations (QSSA) to simplify the dynamics of intermediate complexes such as \(C\) (the CcdA--CcdB complex) and \(C_2\) (the ternary complex). In the above equations, \(A\) and \(B\) represent the concentrations of the antitoxin (CcdA) and toxin (CcdB), respectively, while \(A_{\text{mRNA}}\) and \(B_{\text{mRNA}}\) denote their corresponding mRNA transcripts. The terms \(k_{\text{bind}}\), \(k_{\text{deg}}\), and \(k_{\text{trans}}\) serve as analogues to the Michaelis--Menten constants that control the rates of binding, degradation, and translation. Specifically, the interaction term \(k_{\text{bind}}AB\) represents a bimolecular association that parallels the enzyme--substrate binding step in the canonical Michaelis--Menten model: \[ E + S ⇌[k_{-1}]{k_1} ES \xrightarrow{k_{\text{cat}}} E + P. \] In our system, \(A\) and \(B\) play roles analogous to enzyme and substrate, forming a transient complex \(C\), which may further associate to form \(C_2\). The degradation and regeneration terms (e.g., \(k_{\text{Cdeg}}C\), \(F k_{\text{Adeg}}C_2\)) then capture the turnover and dissociation processes akin to product formation and enzyme recycling in enzyme kinetics. The assumption of a quasi--steady state implies that the formation and breakdown of complexes (\(C\), \(C_2\)) occur much faster than the production and decay of mRNA or free proteins, allowing us to treat the complex concentrations as approximately constant over the timescales of interest. This greatly simplifies the mathematical formulation while preserving biologically realistic dynamics. Because Michaelis--Menten kinetics has been extensively validated across a wide range of biochemical systems, it provides a natural modeling framework for the CcdA/B interaction network. It captures the saturable, nonlinear relationships between toxin and antitoxin concentrations and allows for parameter estimation based on literature data (such as those from the 2020 BITS Pilani iGEM team). In doing so, the model balances mechanistic fidelity with computational tractability, enabling the exploration of steady--state and transient behaviors under various expression regimes.

Python Code

As part of this section, I would like to discuss an important assumption we have made that ultimately governs how our model works. We have modeled the system as a first-passage-time system for the purposes of answering our first question in our reasoning section. We assume that after CcdB expression first overtakes CcdA expression, the cell's growth will arrest, and it will eventually die. This allows us to find the average kill-switch activation time and further extend it to cell death, which we believe to be a reasonably consistent length of time after growth is arrested, of approximately 3 hours, as found by "Gyrase inhibitors induce an oxidative damage cellular death pathway in Escherichia coli" (Dwyer et al., 2007)[4]. Practically, our entire model has been coded in Python and using UVA's High Performance Computer (HPC) Rivanna, and there have been a few very important features that we have made great use of. All of the population dynamics modeling would not have been possible without the use of SLURM, which is an open-source parallel computing software commonly used for node allocation in HPCs. Using SLURM batch scripts, we were able to run our model Thousands of times in parallel and receive large datasets over which we could average. Additionally, for our mechanistic model, we made use of the following open-source Python packages

• Matplotlib
• Scipy
• Math
• Numpy
• Pandas
• Sys

The documentation for which is cited under attributions. For solving our systems of coupled ODEs, we specifically used Scipy integrate with the BDF method because of the stiff nature of our systems. Lastly, all of the code used to generate the model can be found at https://github.com/dsd5kq/UVA-iGEM-2025-Model-

Results

For this section, I would like to recollect the guiding questions that we started with

• How long on average does the killswitch take to activate?
• What is the average internal concentration of cysteine at that time?
• What is the desired promoter of ccdA such that we ultimately reach our target concentration at cell death?

and to answer these questions with the data we have collected. I will start with question 1, which is the main output from the mechanistic model. Due to the stochastic expression of CcdB, output graphs can vary drastically; to remedy this, Figure 1 overlays 1000 trials of our system. The black bars represent the average first passing time and its standard deviation with \(\bar{t}_f=19469.38\pm8937.72\) s. Firstly, the standard deviation is notably quite large, though this is expected as the system is known to and expected to have a certain variance. By definition, we have induced variability into this system, and therefore, this cannot be used as an indicator of model accuracy. Though we can clearly see a clustering of values within that range that begins to thin out significantly beyond. To that end, we find that \(64.5\%\) of our \(t_f\) values fall within the domain \(\{0,\bar{t}_f+\sigma\}\). Meaning we estimate that approximately 2/3 of the population's kill-switch will activate within the first 0.5 hours of a significant increase in cysteine concentration. Figure 2 displays this information in a slightly more cohesive manner and displays a very interesting result. We see that, due to the central limit theorem, we actually recover the Poissonian statistical distribution that we originally put in. This fits quite well, showing us that the decision to use Poissonian statistics was justifiable. An \(R^2\) of 0.8885 is fairly consistent with a good fit; typically, fits closer to 1 are considered better. We still must ask the question of how accurate our average is to answer our second question. To do so, we take the standard error of the mean for our results, which is given by Eq. 35 \[ SEM\equiv\frac{\sigma}{\sqrt{N}} \] Where \(\sigma\) is the standard deviation and N is the number of outcomes in our sample space, for our system, we find that \(SEM=\pm282.64\), meaning our average is accurate within \(2.9\%\) of the actual value. With that degree of accuracy, we can take our average to be a fairly good approximation of the average of our system, and it can therefore be used to answer question 2. By looking at our dFBA outputs at our average time point, we see that the average internal cysteine concentration at the time of kill-switch activation is 2.144 mM. This is of a reasonable magnitude for intracellular cysteine under an over-expression pathway as found in "High Levels of Intracellular Cysteine Promote Oxidative DNA Damage by Driving the Fenton Reaction" (Park and Imlay, 2003)[5]. In this paper, it was found that an intracellular concentration of 1.5 mM is typically reached by a cysteine-overexpressing device. This is indicative that our model may overestimate concentrations in general, which is not unexpected, but is certainly a shortcoming of our model. To find what promoter strength we needed for CcdA we first looked at the available Anderson promoters that we had on hand and their relative strength according to the base promoter. We then used a tool from the Salis lab, specifically from "Automated design of thousands of nonrepetitive parts for engineering stable genetic systems"(Houssain et al., 2020)[6], to calculate the rate constant for the base promoter in units of \(\frac{1}{M*s}\), the equation for which is detailed below \[ k_j=\frac{1}{2}\sum_{i=1}^{2}\frac{\frac{RNA_{i,j'}}{\sum RNA_i}}{\frac{DNA_{i,j'}}{\sum DNA_i}} \] Where i represents the replicate trial, j' represents the j'th concentration amongst all of the listed promoters for both DNA and RNA. These measurements were made for 2 replicates and are then averaged over those 2 replicates to increase the accuracy of the final rate constant. From this Eq. 36, we were able to determine the strength of the base level Anderson promoter and, therefore, the strength of the others that we had on hand.

Promoter	Relative strength to BBa\_J23100	Rate (\(\frac{1}{M \cdot s}\))
BBa\_J23100	1	0.1834
BBa\_J23102	0.86	0.157724
BBa\_J23105	0.24	0.044016
BBa\_J23109	0.04	0.007336

Relative promoter strengths and transcription rates.

\end{table*} For each of our 3 available promoters, we tested the above values for the expression rate constant of CcdA, \(k_{AMRNA}\), to see how it affects our dynamical system. The results of this are in Figures 3 and 4: as can be seen in the case of the BBa\_J23102 Figure 3, the level of expression was such that it would always outweigh the expression of CcdB. This would therefore make activation of our kill switch almost impossible, meaning we would have no control over cell death. Additionally, looking at the case of BBa\_J23109, Figure 4, we see that expression of CcdA cannot outweigh that of CcdB. This does not allow sufficient time for cysteine to overexpress and is therefore not viable. With BBa\_J23105, Figure 1, we are in the goldilocks zone, and expression of the CcdA/CcdB systems perfectly outweigh each other, allowing for fine-tuned control of the kill-switch. Ideally, we would prefer a slightly less powerful promoter, as currently, we would reach an external cysteine concentration of about 433 mg/L with the BBa\_J23105 promoter. To reach our desired external concentration of 300 mg/L, we would need an average first passage time of 15945 s. An additional constraint on this is the availability of discretized promoter strengths, necessitating an approximation based on the closest available strength. Based on rudimentary analysis, we have found that the closest existing Anderson promoter to our desired concentrations is BBa\_J23105, but we will discuss achieving the value of the ideal Anderson promoter in Section 5.

Conclusions

From our model, we were able to answer our 3 main guiding questions, which are listed below for accessibility

• How long on average does the killswitch take to activate?
• What is the average internal concentration of cysteine at that time?
• What is the desired promoter of ccdA such that we ultimately reach our target concentration at cell death?

We found that using the best approximation of our desired systems with existing constituitive promoters, the total time to reach kill-switch activation on average is approximately 0.5 hours after cysteine production spikes significantly. Additionally, it was found that the average internal concentration of cysteine at that time would be 433 mg/L, 133 mg/L above our desired target concentration, though this could very well be an overestimate. Given the existing sequences, the best promoter we currently have for the expression of CcdA in our system is the Anderson promoter BBa\_J23105. This existing data provides an excellent testing framework for wetlab data and has many places in which it could easily be informed by wetlab data to improve the accuracy of our model. The model serves as a fairly accurate model of CcdA/CcdB system that could easily be appropriated and reused for iGEM teams in the future, as the code is well commented. Additionally, it provides a framework of modified Hill equation ODEs that could again be used for other systems that are not fully characterized, such as the CcdR system. Lastly, we introduce a somewhat novel approach for modeling the stochastic expression of genes using a 2-state Markov chain and binding/unbinding probabilities based on Poissonian statistics. A method excellent for characterizing the fluctuations that arise in low particle density situations such as ours. This code is also well commented and could easily be adapted for use in other systems, serving as a novel contribution to future iGEM teams in characterizing similar systems.

Future Work

Our future works include many important aspects of the project that could be expanded upon and others that would also benefit the project, but have not yet begun. The first thing we would like to do is to create a curve based on the existing Anderson promoters that would allow us to directly characterize our desired promoter strength based on the final external cysteine concentration. This would allow us to potentially look for other constitutive promoters of the same strength, getting us closer to our desired output value. Looking beyond, incorporating post-killswitch-activation expression of cysteine would also be invaluable, as it would better allow us to characterize how the final cysteine concentration will look. Ultimately, incorporating that into our population dynamics in a system that feeds data from the dFBA to the mechanistic model and back to the dFBA to inform the next step would be the final goal in most accurately characterizing our given system.