Mathematical Modeling of a Bioluminescent E . Coli Based Biosensor

Abstract. In this work we present a mathematical model for the biorepor ter activity of an E. coli based bioluminescent bioreporter. This bioreporter is bas ed on a genetically modified E. coli which harbors the recA promoter, a member of the bacterial SO S response, fused to the bacterial luminescence ( lux) genes. This bioreporter responds to the presence of DNA damaging agents such as heavy metals, H 2O2 and Nalidixic Acid (NA) that activate the SOS response. In our mathematica l model we implemented basic physiological mechanisms such as: the penetration of the NA into the biosensor; gyrase enzyme inhibition by the NA; gyrase level regulation ; creation of chromosomal DNA damage; DNA repair and release of ssDNA into the cytoplas m; SOS induction and chromosomal DNA repair; activation of lux genes by the fused recApromoter carried on a plasmidal DNA; transcription and translation of the lumin escence responsible enzymes; luminescence cycle; energy molecules level regulation and the regulation of the O2 consumption. The mathematical model was defined using a set of ordinary dif ferential equations (ODE) and solved numerically. We simulated the sy stem for different concentrations of NA in water for specific biosensors concentra tion, and under limited O2 conditions. The simulated results were compared to experim ental data and satisfactory matching was obtained. This manuscript presents a proof of c oncept showing that real biosensors can be modeled and simulated. This sets the g round to the next stage of implementing a comprehensive physiological model using experimentally extracted parameters. Following the completion of the next stage, it w ill be possible to construct a “Computer Aided Design” tool for the simulation of the gene tically engineered biosensors. We define a term “bioCAD” for a Biological System Computer Aided Design. The specific bioCAD that is described here is aimed towards wh ole cell biosensors which are under investigation today for functional sensing. Usag e of the bioCAD will improve the biosensors design process and boost their performance. It will also reduce Non Recurring Engineering (NRE) cost and time. Finally, using a parameterized solution will allow fair and quick evaluation of whole cell biosensor for various applications.


Introduction
Today, the design of biosensors does not include a computer based simulation stage as a standard procedure similar to what was a standard use in most engineering disciplines.The current commonly-practiced "trial and error" method which is used for biosensor design is a complex, expensive and time-consuming process that fails to guarantee the required system performance.Therefore, currently there is a need for a computer aided design (bioCAD) tool for biosensor engineering.
Any CAD tool is based on mathematical modeling of the resources and processes of the corresponding object.The bioCAD for the whole cell biosensors engineering has to cope with cellular functions that are based on a large number of chemical reactions and transport processes that are regulated by proteins, enzymes, genetic and biomechanical mechanisms [1][2][3][4].The relatively simple mathematical models deal with enzymes and substrates by ordinary differential equations (ODE) [5][6][7][8].A more sophisticated form, considering diffusion, is based on the partial differential equation (PDE), as shown for example in works describing a biosensor acting in a trigger mode [9] and in the Vibrio fisheri control system [10].The genetic and biomechanical mechanisms could be found in works depicting the E. coli SOS function [11][12][13] for damaged DNA repair.The modeling could be stochastic or deterministic, depending on a number of substrate and reactant molecules [14][15][16] that are being taken into consideration.
Metabolism is commonly defined as a process in which nutrients are converted to provide energy and the synthesis of new organic materials for the cells maintenance activity and reproduction [17,18].Therefore, metabolism modeling usually appears in most whole cell biosensors [16].This modeling is usually very complex and it is common to derive algorithms for simplifying the metabolic networks [6,19,20].The second major issue is the modeling of the regulatory networks.For example, those were comprehensively studied and modeled for E. coli [6,16,19,21], and for Synechocystis [22].
The other main parts of the biosensor mathematical model are: (i) the interaction of the analyte with the biosensor; (ii) the biosensor response; and (iii) the reporting mechanisms.The analyte interaction with the biosensor, which is the initial link in the detection chain of the bioreporter, has been reviewed in the literature, for example see references [23,24].The biosensor response activates the reporting mechanism.The SOS mechanism (a response tool for the repair of a DNA damaged by the analyte) is described in a variety of works [11-13, 25, 26].Other response mechanisms models are based on transcription of signaling proteins and enzymes [27,28].In our research we are investigating luminescent bioreporters that are based on the activation of bacterial luminescence genes (lux) derived from marine luminescent bacteria Vibrio fisheri, all of them are well-studied mechanisms [29][30][31][32].
We use "normally-off" E. coli based bioluminescence biosensors incorporated into 1-100 µL reaction chambers, densely packed on a biochip.Special µFluidics on the biochip are used for translocation of the tested liquid with different concentrations of Nalidixic Acid (NA) into the reaction chambers.O 2 supply, which is vital for bacterial metabolism and for bioluminescence, is relatively limited.The absence of excess O 2 supply makes the biosensors environment "unfriendly" and limiting the light generation efficiency.
The implementation of Stage-I of the biosensor Computer Aided Design (bioCAD) roadmap, see Fig. 1, is demonstrated in this work.The roadmap consists of three stages.The Stage-I is for bioCAD proof of concept for the whole cell E. coli based bioluminescent biosensor.Only the main physiological processes are used in Stage-I.In this stage most of the reaction and process parameters are taken from the literature or, alternatively, assumed by some intelligent guess which is appropriate to the physiological mechanism.The simulated biosensor behavior matches in general the kinetics from the experimental results and therefore lays the ground to Stage-II.During Stage-II of the roadmap (Fig. 1), the mathematical model should be expanded to a "Comprehensive Physiological Model" while the parameters should be experimentally measured for the specific biosensor used in the research.For the successful completion of this stage the simulation should match exactly the experimental behavior.Once completed, the model could be used for the computer aided design to improve the biosensor performance which is the main part of Stage-III (see Fig. 1).

Reduced physiological model
There is a diversity of factors that could be desired to boost the performance of the biosensor.Here a few are mentioned: (i) shortening the response time between the analytes incursion and the reporting initiation; (ii) increasing the reporting intensity; (iii) defining a selective response for specific analytes; (iv) adjusting the amplitude of the reporting response corresponding to the analytes concentration, etc. Design and simulation of the new biosensor will tremendously reduce the demands for time, labor and financial resources.Once the desired performance is met, the actual biosensor can be built and applied.Next, the bio sensor's experimental behavior should be compared to that of the simulation, and the difference between them should be studied.A good model will yield a reasonable match, yielding a working functional biosensor applicable for system interfacing.Under other circumstances, it is necessary: (a) to check the physical implementation of the biosensor and (b) to check the comprehensive physiological model and repeat the evaluation process until the matching is achieved.
In this section we have presented a diversity of work dealing with modeling of various biological processes.However, the complete physiological model of genetically engineered biosensor based on E. coli, SOS response and luminescence reporting has not been implemented yet.In the current work we did such an implementation, successfully validated it, and now it is possible to proceed with the 2nd step toward biosensor computer aided design development.

Biosensor mathematical model
In this work we build the complete model of a genetically engineered bioluminescent E. coli based biosensor using subsets of equations describing the critical physiological mechanisms.Note that, as described before, only to the conceptual model is referred to, where only the rate defining steps are taken into consideration.Simplified models for those steps have been used.
Those mechanisms are (i) analyte uptake and SOS response; (ii) simplified metabolism regulation model; (iii) simplified O 2 consumption; and (iv) light generation process.All the aforementioned mechanisms have mutual dependency creating a network or a biological feedback circuit.Activation of each mechanism depends also on the physiological state of the biosensor.The physiological states of the biosensor are described in the following sub-section.

Biosensor's physiological states
The biosensor under investigation referrers to bacteria that are positioned within reactions chambers on a chip.In this case the biosensors are typically enclosed within porous polymer matrix, such as oxide Sol-Gel or agar, and we may assume that mitosis is prohibited.Therefore, the possible main physiological states are as shown in Fig. 2.  Initially the biosensor is in an "Idle State".Once analytes are introduced, the internal machinery of the biosensor is damaged and it becomes under a "Stress" condition.Following the stress, the sensor enters into a "Recovery State".In our particular case the stress is actually DNA damage and hence the "Recovery State" is combined of two processes: (i) a SOS response and (ii) a "Metabolic Recovery".If the stress induced injury is above some threshold level the cell dies, i.e. enters the "Cell Death State".Otherwise it generates a reporting response, i.e. enters the "Response State".The processes of stress, metabolic recovery and SOS, could continue while the sensor is in the "Response State".Once the stress and the recovery processes are completed the sensor returns to the "Idle State".

Analytes uptake dynamics and SOS response activation
In this section the process of analyte uptake by the biosensors and the activation of the SOS response are shown.First the problem is reduced to a single biosensor.The model simplification is achieved by assuming uniform distribution of the biosensors and analytes inside the reaction chamber.Each biosensor is treated as a "reaction center" as demonstrated in Fig. 3. Assuming a reaction chamber with volume V RC and the biosensor concentration of C b , the number of the biosensors inside the reaction chamber N b is: We also assume that the average volume of reaction chamber for single biosensor, V RC/b is: In this case, the average distance, L b2b , between two adjustment biosensors is: We assume that the concentration of the introduced analyte is uniform through the reaction chamber and each biosensor, in average, shares the same amount of the analyte molecules.The analytes concentration, C an−ppm , is given in "parts per million" (ppm) that is equivalent to mg/L.The average number of analyte molecules N an/b per partition sub-volume, V RC/b is given by: where M W an is a molecular weight of the analyte and N A is an Avogadro number.For example, for 1 ppm of Nalidixic Acid (NA) with M W NA = 232 gr/mole, reaction chamber with volume V RC = 1 µL, biosensors concentration C b = 5e5 b/µL, the number of the NA molecules per volume of reaction chamber utilized by a single bacterium is: The Nalidixic Acid (NA) analyte emulates a group of toxicants that damages the chromosomal DNA of the bacterium.Note that NA is relatively safe for humans and, therefore, preferably used to test the performance of such biosensors; consequently we selected using NA in our work.The simplified scheme of the NA uptake and SOS process activation is depicted in Fig. 4. We assume that the analyte transport in the liquid (free or enclosed in a porous matrix), via diffusion or convection, is not the rate determining factor.For example, taking the data in the abovementioned test case, the average dimension of the reaction chamber volume utilized by the biosensor is The diffusion coefficients of most water-soluble compounds with low molecular mass are in the order of D ∼ 1e − 5 cm 2 s −1 [24].Therefore, the average diffusion time of such molecules, t b2b is The uptake rate of the analyte like NA is faster in few orders of magnitudes than the resulted reaction (0.15 sec versus minutes).Therefore, the transport of analytes like NA is not limited by diffusion and it is not taken into consideration.The analytes penetrate into the bacterium through uptake channels available on the membranes of the biosensor.The change of NA concentration inside the biosensor is expressed by the following equation: where k trn−NA is a concentration gradient coefficient indicating the percentage of molecules transported through the cell membrane per second, NA in , NA out are the NA concentration inside and outside the biosensor correspondingly, H NA is the half life time of the NA molecule degradation, k GyrA.NA is a constant rate of interaction between NA and gyrase enzymes whose concentration is marked by GyrA.We also use molar concentration, NA inMolar that is converted from the number of molecules per bacteria using the Avogadro number N A and the bacterium volume V b : Outside the biosensor the change in NA concentration is as following: The NA analyte inhibits the gyrase enzymes.The gyrase enzyme is responsible for unwinding supercoiled DNA during replication.However, when gyrase is attached to the DNA, but inhibited by NA, it fails to function and DNA lesion is created in that region [33].In our model we assume that DNA polymerase enzymes repair the dimers.The repair is modeled by polymerase motion on the chromosome which causes the removal of the lesions.We will relate to the lesion's site as a "DNA dimer".Once the lesion is detached and corresponding DNA fragment is repaired, the dimer is removed.The change in the number of dimers, N dim , could be expressed as: where dN dim−gen and dN dim−rep are the number of dimers generated and repaired correspondingly.The number of the generated dimers depends on the number of inhibited gyrase enzymes GyrA inh attached to the chromosomal DNA (GyrA inhDNA ).The rate of damage production is expressed by means of a rate coefficient k dim : The concentration of the inhibited gyrase molecules attached to DNA (GyrA inhDNA ) is mediated by a constant k GyrA.DNA , the gyrase half life time H GyrA and also decreases through the dimers repair process: The change in the total number of the inhibited gyrase enzymes is a function of the NA ⇔ gyrase binding coefficient k GyrA.NA : The cell can regulate gyrase concentrations [34].In our model we assume that the gyrase concentration decreases due to NA inhibition and some effect on transcription.The gyrase concentration and production rate deviates from their normal values, GyrA norm , and T GyrAnorm , respectively.In this case the rate of gyrase production, T GyrA is accelerated.This rate is modeled by the following function of some arbitrary coefficients GyrA min and P GyrA , and also GyrA norm which was described before.
Next, we present the change in the gyrase concentration that depends on the gyrase production rate T GyrA , the inhibition of the gyrase by the NA and the half life time The dimer removal rate depends on the average distance (in nucleotides) between the damaged DNA, L d2d , the velocity of DNA polymerase (DNAP) movement on the DNA, v s , and the delay, t dim , that takes to repair the damaged part of the DNA: where the DNAP is the number of DNA polymerase enzymes in a cell and k bind is a binding probability of those enzymes to the chromosomal DNA.The coefficient 2 that appears in equation ( 13) is due to the two replication forks that are moving in opposite directions at similar rates.The average distance between two dimers is a ratio between the lengths in nucleotides of chromosomal DNA (L ChrDNA ) to the current number of dimers: The DNA repair process or SOS response detaches the damaged DNA parts.As a result, ssDNA fragments are released to cytoplasm.The RecA molecules that are normally under stable, low concentration bind the ssDNA fragments and are transformed into protease enzymes with constant rate k protease .The protease enzymes interacting with LexA proteins cause the autodigestion of LexA with a constant rate k LexA .The LexA protein is a direct repressor of the lexA and recA genes therefore the rate of synthesis of the LexA and RecA proteins is proportional to the DNA template of the lexA/recA locus that is free from the repressor [11].In the absence of the repressor, the LexA production rate is T LexA and the RecA production rate is T RecA , while EK LexA and EK RecA are the corresponding equilibrium constants.The changes in the concentrations of the ssDNA, P rotease, LexA and RecA, considering the H ssDNA , H P rotease , H LexA and H RecA half life time constants, are expressed by a set of the following equations: We just have shown the SOS process where the kinetics of the LexA protein concentration depends on the Nalidixic Acid analyte inducer concentration (C an−ppm ).Later we show how the LexA level influences the intensity of the bioluminescence.

Metabolism as a function of DNA damage
We found, experimentally, that upon exposure to NA concentration above a certain level, the initial (background) luminescence activity decreases.Our assumption is that the DNA damage, that initiates the SOS process, also disturbs the rate of normal-baseline metabolism.This baseline metabolism is also responsible for the generation of reducing power and energy molecules like ATP, FMNH2 and NADPH, which we will refer to as EN RG, (Fig. 5).In our model we express the ENRG concentration as a function of the number of dimers N dim using the following arbitrary fitting function:

Metabolism Enzymes
where, P dim is a fitting coefficient.The ATP, FMNH2 and NADPH are also necessary for the light generation process.The ATP is also a driving force for the DNA polymerase.Therefore, the decrease in ENRG concentration is also reflected in the velocity of DNA polymerase movement on the DNA also expressed using arbitrary fitting function: where P EN RG−v is a fitting coefficient and v s−max is the maximal DNA polymerase velocity.The DNA polymerase promotion velocity influences the rate of DNA repair and ssDNA generation, as described in the SOS section.

Molecular oxygen consumption
The biosensors described in the current work are encapsulated into the reaction chambers of the biochip [35].In the current system configuration, the biochip is sealed and therefore a constant oxygen supply is not obtainable.The only the available oxygen for the reaction is present in the sample to be analyzed.Normal concentration of dissolved oxygen in water is CO 2 given in [mg/L] units.It corresponds to molar concentration of: where, M W O2 is a molecular weight of the molecular oxygen.The number of the O 2 molecules per liter, N O2−ltr is: and the number of the O 2 molecules per bacterium is: We assume that O 2 freely diffuse through the bacterial membranes inside and therefore the internal O 2 concentration is equal to that of the medium of the reaction chamber.The oxygen consumption by bacterium in an aerobic condition is K O2−mol , [36] where it is expressed in [moles/min/cell] units.We transform this into [molecules/sec/cell] units and get K O2norm : The O 2 is consumed by the bacterium in two pathways.The first is due to normal aerobic metabolic activity and the second is due to the bioluminescence reaction.We assume that bacteria can regulate its level of aerobic-anaerobic activity and decrease the level of O 2 consumption correspondingly to the O 2 concentration fitted by P O2uptake coefficient.Therefore, the change in O 2 molecules available per bacterium could be expressed as follows: where N photons is the number of photons emitted per dt time and QE luc is a quantum efficiency of the bioluminescence process.

Bioluminescence generation process
The bioluminescence generation process implements the reporting part of the biosensor.The detailed description of the process presented in Appendix, while in this section we implement a simplified model, see Fig. 6(b).Bioluminescence generation occurs due to the presence of enzymes generated by lux genes and a number of substrates participating in a "Light Generation Cycle" (LGC), as shown in Fig. 6(a).The substrates required for the LGC are: (i) energy molecules: AT P , NADP H and F M N H 2 along with (ii) other molecules: O 2 , H 2 O and activated fatty acyl groups (RCO.ACP ).In this work we assume that substrates are in excess, except for energy molecules and O 2 .
There are two groups of the lux related enzymes participating in LGC.The first group is a complex of three enzymes: transferase, synthetase and reductase, presented as triangular in Fig. 6 with letters "t" "s" and "r" correspondingly.The primary function of those enzymes is the generation of luciferin from the substrates.The second group consists of single enzyme luciferase, "l", responsible for the generation of the light emission from the luciferin.
The induction mechanism, indicated by an arrow "Promotion" in Fig. 6(a), is different for each kind of the promoter.For the recA promoter, used in the current work, the promoter activation is repressed by the LexA protein bound to the recA activation site.Therefore, the promotion intensity is a function of the LexA protein concentration in the cell, described in the SOS sub-section.We represent the promotion intensity, I P r , by the following expression: where I P r−min is the basic promotion intensity and LexA norm is a concentration of LexA for the non-induced state of the biosensor.The LexA min is an experimental parameter taken from the literature [11,37].The enzymes are generated through plasmidal DNA transcription and then the translation process.The number of plasmids in the cell, N plasmids , is controllable during the biosensor genetical engineering.First, mRNA containing luxABCDE coding sequence is transcribed with the maximal rate of T c , multiplied by the number of plasmid copies and by the promotion intensity factor.The mRNA concentration, considering the mRNA half-life time H RNA , is expressed by the following differential equation: Next, the enzymes are translated by ribosomes moving on the mRNA patterns with an enzymes generation rate of T l .Enzymes, like all protein molecules, are going through degradation processes and are therefore characterized by a decay rate or half-life time.The luciferase enzyme is considered as a stable protein with half life time H sp , while the transferase, synthetase and reductase are unstable proteins with half life time of H up [38].The concentrations of the luciferase E l and reductase, synthetase, transferase E r,s,t are expressed by the following differential equations: The simplified LGC, see Fig. 6(b), consists of enzymes, external inputs, intermediate products and output as photon emission.The intermediate product Fatty Acid (or RCOOH), whose concentration is signed by F A, is created by transferase E t and luciferase E l enzymes activity.The Fatty Acid is converted to Luciferine, whose concentration is marked by L, by means of synthetase and reductase enzymes E sr .The conversion process requires also ATP and NADPH energy molecules, whose concentration is signed by EN RG.Finally, luciferase converts the Luciferine molecule into photon.This reaction requires also molecular oxygen O 2 and F M N H 2 energy molecule.
The change in the concentration of Fatty Acid depends on the reaction constants that relate to transferase, k t , luciferase, k l and synthetase-reductase, k rs , activities and Fatty Acid half life time H F A : The luciferine is generated from the fatty acid in a chain reaction mediated by synthetase and reductase enzymes and also depends on consumption by luciferase and the luciferine half life time, H lcf r : The k t reaction coefficient relies on reaction constant k t−basic and also treats the situation where the product is above saturation concentration F A sat : The generated product fatty acid is used as a substrate for luciferine product generation mediated by synthetase-reductase, k rs , reaction coefficient.This coefficient is expressed through reaction constant k rs−basic , adjusted by substrate saturation concentration F A sat , and availability of the energy molecules: The same is for luciferase related reaction coefficient k l , while it is also adjusted by availability of the molecular oxygen.
The number of the generated photons per unit time is proportional to the k l constant and number of the luciferine molecules along with luceferase enzyme and quantum efficiency of the light emission QE luc .The number of the emitted photons N ph is expressed in the next manner: We have expressed a photons emission process that depends on LexA, EN RG and O 2 input levels.This section completes the puzzle of modeling the physiological processes of the genetically engineered biosensor described in the current work.Following, it is shown how this model is implemented as a computational program.

Numerical implementation of the mathematical model
In this section we show an implementation of the biosensor response mathematical model.First, we show the general technique of the numerical solution using iterations, following reaction constants and initial conditions are presented, and finally the outcomes of the simulations are shown.

Numerical solution by iterations
We implemented the modules of the mathematical model, described in the previous section, using numerical method.First, we defined interaction step dt and a number of iterations according to the total simulation time requirement T sim , that is usually no longer than 5 hours.Next, we ran loops? of n intr iterations: Most equations used in the model are expressed in the form of Ordinary Differential Equation (ODE): dX(t)/ dt = F (X(t)).In order to implement them numerically, we calculate the change dX for single iteration i, and the value of the X for next iteration - We have to use max(0, X(i)+ dX) because X represents a concentration or number of molecules that cannot be negative.The precision of the equations solution, as well as run time, is a function of dt value.We found that dt of 1 sec provides a smooth enough solution and reasonable run-time (about 1, 000 iterations per minute on 1.6 GHz PM processor).

Experimental inputs, initial conditions, reaction coefficients and variables
For conversion of the mathematical model developed in the previous section we need to substitute initial conditions and coefficients for the equations presented there.Parts of the values are available from references.The other part was not found in references and either should be experimentally measured or assumed.In this work missing data was estimated based on expected values and adjusted to yield a good behavior.
In this subsection the values used in the mathematical models are gathered into the tables.First, the experimental inputs are presented in Table 1.Next, intermediate variables are gathered into Table 2. Then we put initial conditions used for the variables into Table 3.After that, we show the fitting coefficients, see Table 4, and finally, rates of synthesis, reaction constants and coefficients are placed into Table 5.

Simulation results
In this section we present results of simulation runs of the mathematical model implemented in this work.We used reaction coefficients and initial conditions mapped in the previous section.The simulation input was set according to the experiment presented in the next section and it corresponds to the typical values in Table 1.The following concentrations of the NA analyte were used: 0 ppm, 0.78 ppm, 3.13 ppm, 12.5 ppm, 25 ppm, 50 ppm and 100 ppm.At the first stage the simulation takes the biosensors to the steady state.This steady state emulates preparation of the biosensors under ideal incubating conditions without limit in O 2 , nutrients and absence of toxic materials in the medium.The process lasts about 150 minutes.The luminescence reaches a stable base level making it possible to pass to the second stage of simulation, where the biosensors are encapsulated into the observation chambers without additional O 2 and nutrients supply and are challenged with various concentrations of the NA analyte.
The simulation results for the second stage are presented in Figs. 7, 8.Each figure contains sub-plots exhibiting kinetics of components used in the model.The components are identified according to the title on the top of the sub-plot.The first, Fig. 7, shows results of the SOS process; while the second, Fig. 8, demonstrates kinetics of the light generation cycle.The x-axis for each sub-plot is time in minutes.In the next section actual results for the same Nalidixic Acid concentrations, as in simulations, are presented and compared.The experimental results contain only kinetics of emission intensity versus time.Therefore it is extremely challenging to analyze bottlenecks for experimental results, while, as could be seen from this section, using simulation is very friendly.

Experimental results and discussions 4.1 The experimental setup
The measurements were performed on a Victor-II luminometer by Wallac Inc. using a 384 chambers microtiter plate.To avoid cross-lighting, each second chamber was utilized for the experiment.The working chambers were fully filled by a solution containing Luria Bertani (LB) medium, bacteria under concentration of 5e8 bacteria/mL and different concentrations of Nalidixic Acid (NA).The working chamber volume is about 50 µL, and therefore included about 2.5e7 bacteria.The top of the microtirtle plate was sealed by a transparent film to avoid oxygen supply (in order to emulate biochip condition).We used the following NA concentrations: 0 ppm, 0.1, 3.13, 12.5, 25, 50, 100 ppm.In addition, chambers with water only and with water and LB only were prepared.Four repetitions were set up for every NA concentration.The light intensity from the each reaction chamber was sampled once in 3 minutes.

Experimental results
The experimental results are demonstrated in Fig. 9.The y-axis is the total collected light intensity in [photons/sec] units and the x-axis is the time since analytes intervention in minutes.There is an individual curve for each NA concentration.The experiment was performed during 240 minutes similarly to the time of the simulation.

Discussions
As can be seen, the kinetics of simulation matches quite reasonably the experimental results (see Fig. 10 and Fig. 9).The biosensors are in Stress, Metabolic Recovery and SOS states (see Section 2.1) for the first 45 minutes.The luminescence intensity reaches the same order of magnitude 1e7 [photons/sec] during the "Response State" and finally the luminescence activity discards after about 180 minutes correspondingly to the reversion of the biosensor to the "Idle State".Even though the kinetics as a function of NA concentration did not match precisely, the main goal of the current work, Stage-I, has been achieved (see Fig. 1 in Section 1).We still have to note that the quality of results during the Stage-II can't be fully estimated.This is due to complexity of the presented model from one side and simplifications of the biological processes from the other, as well as possibility of non-unique solutions for the fitting of experimental data.

Conclusions
In this work we have presented a design of the mathematical model of the response of a whole cell based biosensor, as well as simulation results generally matching the simulation.We have defined the three stages for implementation of Biological Computer Aided Design (bio-CAD) software and successfully fulfilled the first stage.The completion of the next two stages would require collaborative work of (i) bio-mathematical modeling; (ii) molecular biologists; and (iii) engineers.This extremely challenging task, considering obstacles such as complexity of the model and non-unique solutions for the inverse problems and may demand 10 to 20 man years of R&D.However, once the bio-CAD software would be available, the development of the whole cell based biosensors will be boosted in a few orders of magnitude.

Appendix. Comprehensive model of light generation cycle
In this appendix the information from the prior-works [5,30,31,42] is integrated into a comprehensive LGC chemical and mathematical model.This process is illustrated in the chart below.The lux enzymes E l , E r , E t and E s are generated by corresponding genes luxAB, C, D, E with the generation rates k AB , k C , k D and k E .Those generation rates are different to each gene and have units of [M/sec/Pr].The enzymes are marked by a triangular shape.The enzymes are also having half life time constants k e0.5τ .The transaction between the reactions compounds are accompanied with the chemical reaction rate constant k rate (see Table below).Part of the reactions are reversible and their constants are marked k −rate correspondingly.

Fig. 1 .
Fig. 1.Stages for development of computer aided design for bioluminescent E.Coli based biosensor.

Fig. 4 .
Fig. 4. NA uptake and chain reactions of the SOS process initialization.

s
The names of molecules participating in the LGC are presented in the

Table 3 .
Initial conditions

Table 5 .
Rates, constants and coefficient