Team:Imperial College/Model
Modelling
Diffusion modelling
Analytical Evaluation
We wanted to perform the initial evaluation for the concentration profile by hand to ensure we knew every assumption made in its entirety, and to maximise our understanding of the system. We did not want to simply use numerical modelling software without first understanding the key concepts and limitations within the system to be able to identify aberrant results and unjustified assumptions.
Figure 1. Illustration of the PixCell Construct
We initially modelled the diffusion of the chemical species away from the electrode by treating the electrode as a disc source with a constant flux of the molecule across it. We then constructed the problem by assuming that the agar medium itself was not generating additional molecules, and therefore that the only source of said molecules was the disc electrode. This assumption, combined with the conservation of species, allows us to describe the concentration of the species within the agar with the partial differential equation known as the Laplace equation. The Laplace equation in this situation represents the idea that the number of molecules arriving at a point in space per unit time is equal to the number of molecules that leave that space per unit time. The radially symmetric nature of the problem (the problem looks the same from any angle, as the electrode is a disc, not a rectangle or any other shape) allows us to use the 2D axisymmetric form of the Laplace equation, rather than the full 3D form, which makes solving the problem easier.
Due to the slow speed of diffusion, and the square root proportionality between time and diffusion distance, we assumed that the agar plate was arbitrarily deep, and that the concentration of the species tended to 0 as depth (z) increased. This can also be conceptualised as the diffusion “front” not reaching the bottom of the plate within the timescales that we are interested in. This assumption was also applied in the radial direction (ie the diffusion front does not reach the edge of the agar plate). This assumption was later confirmed through the use of Comsol modelling. (INSERT LINK?) We also constrained the problem to positive values of z, meaning that no molecules can leave the agar through its top surface. We also assumed that the chemical species analysed would not adsorb onto, or otherwise interact with, the agar medium.
The full mathematical solution for the problem, as well as the justification for the assumption that the molecules are only transported by diffusion, is given in the appendix However, we shall say here that the problem was formulated in cylindrical coordinates (again due to the radially symmetric nature of the problem), and non-dimensionalised both to simplify working and to generalise the problem to any desired scales (scales here meaning length scale of the problem and scale of diffusion, normalised for diffusion constant and flux through the electrode). The scalings are given below. The problem was evaluated analytically and by hand up until the last step. That last step involved the evaluation of a complicated infinite integral, which was non-trivial. However, because the integrand tends to 0 as the variable tends to infinity, use of numerical integration in this case was justified to give the solution for the dimensionless concentration profile. SCALING FORMULAE INSERT (IntegrandPlot1) The exact magnitude of concentration for a given species at a specific flux can be found through simply reinserting the diffusion coefficient for that species. Interestingly, the diffusion coefficients for the species of interest are all reasonably similar, and so give reasonably similar results. Despite this, a wider range of diffusion coefficients was tested (between 5e-10m^2s^-1 and 1e-9m^2s^-1) to account for any differences from literature value due to interactions with the agar medium, or differences due to molecular interactions within the solution.
Observing the scalings for the non-dimensionalised problem, it can be seen that the width of the concentration distribution is directly proportional to the electrode radius, and therefore scales linearly with it. The concentration at any given point is directly proportional to the flux through the electrode, and inversely proportional to the Diffusion coefficient. The half-maximal width of the concentration profile is 1.145 times the electrode diameter, which confirms the feasibility of generating spatially resolved regions of oxidising or reducing environment, with the resolution limited only by the electrode radius, which can be considered an engineering parameter. INSERT JUSTIFICATION FOR SPACING PLOT OF DIFFERENT D (SSforDs) SURF OF PROFILE RADIAL PLOT
Comsol Multiphysics
Comsol is a very high-level Computer-aided design (CAD) tool, which supports the analysis of many types of problems with wide-ranging physics, and the coupling of different physics interfaces together to compute intricate and complicated problems. It has built-in geometry and meshing tools, which allow the user to create model geometry in-package, without the need of external CAD software, and was therefore perfect for performing our simulations. The manual calculation of the diffusion profile, although useful, is confined to steady state, and cannot provide the same level of tunability as a model developed in a dedicated simulation software package such as Comsol. Building a model in Comsol allowed us to check the results of the analytical calculations, to observe the evolution of the concentration profile over time, and to quickly and easily analyse different geometries and test the assumptions made when formulating the analytical problem.
DO WE NEED A HOW? MAYBE A QUICK COMSOL GUIDE??
This model was constructed in order to test assumptions used in the analytical solution, specifically the assumptions that the diffusion “front” would not reach the edges of the plate within the relevant timescale. To do this, we used a time-dependent study with the 2D axisymmetric geometry pictured below, where the disc on the top surface is the source electrode. INSERT SIMPLE 2D AXISYMMETRIC FIG Results and Conclusions As can be seen from the figure below, after 20 minutes, the concentration front was nowhere near the bottom of the plate, let alone the sides of the plate. Note that the edge of the domain shown is not the edge of the plate; the domain shown has a radius of 1cm, and the true edge of the plate would be much further from the electrode than this. The model corroborates the analytical solution, showing a very steep concentration gradient within (INSERT)mm of the edge of the electrode, again making a good argument for the achievement of spatial resolution in a patterning application such as ours. INSERT 20MINS CONC PROFILE (2D SIMPLE) From this model, we can conclude that the assumptions in the analytical solution regarding the diffusion “front” are valid, and that it is valid to treat the dish in the analytical solution as both arbitrarily deep and wide. Additionally, it is reasonable to think that localised oxidised or reduced volumes can be induced in the agar plate in a period of 20 minutes by the use of electrodes, which was confirmed experimentally.
When a manufacturer produces a PCB with pads on the underside, such as the one we ordered, the pads are plated, and then the underside of the PCB lacquered. This means that the actual surface of the pad is slightly recessed from the surface of the rest of the PCB. INSERT IMAGE OF RECESSED BOIIIIIIIIIIIIIIIIIIII EXTRA RECESS LIKE THE GREAT RECESSION Therefore, we wanted to test how a small recession such as this would affect the concentration profile around the electrode. We again used 2D axisymmetric geometry, shown below, where the electrode disc is the top surface of the recessed cylinder. Results and Conclusions INSERT GRAPHIC The simulation shows that even very small recessions drastically reduce the diffusive flux away from the surface. This is because when the surface is recessed, even by a few um, all flux of the species must be perpendicular to the surface. This constraint drastically limits the space into which the molecules can diffuse, as there is a concentration gradient only in the z-direction, as opposed to in both the r- and z-directions in the flush case. This in turn limits the diffusive flux into the surrounding environment, and results in an increased concentration at the electrode surface due to limited diffusion out of the recession. INSERT GRAPHIC- FLUX LINE COMPARISONS It should be noted that the decreased diffusive flux into the wider environment is not necessarily a bad thing, as although it limits the volume a particular electrode can affect, it increases the sharpness of the concentration profile (and hence gives a sharper on/off transition for induced cells at the boundary), and increases the concentration directly underneath the electrode, resulting in greater fold induction of any electro-inducible organisms in this region. Indeed, the magnitude of recession of electrodes within their mounting could be approached as an engineering specification, with a smaller recession resulting in wider but less well-defined fields of influence, and a larger recession resulting in smaller, sharper and more highly induced fields of influence.
Transcription-Translation Modelling
Comsol Multiphysics
In order to better understand the system we are engineering, we constructed a deterministic transcription-translation model of the gene circuit using ordinary differential equations in Matlab. This model allowed us to explore the dynamics of the system, to better understand the parameters that the model would be sensitive to, to use curve fitting tools to estimate unknown parameters, and to check whether the constructs were behaving as expected.
Using parameters from literature, we constructed a model consisting of 4 coupled ODEs: one for transcription of soxR, one for translation of soxR, one for transcription of GFP, and one for translation of GFP. SoxR production was modelled as constitutive, as there was insufficient evidence in literature to support feedback control of the transcription of soxR. The binding of soxR to its promoter was modelled using the phenomenological transfer function given below. There was no clear data to be found regarding the cooperativity of the soxR transcription factor; however, it is known that soxR forms dimers, so the maximum cooperativity is constrained at 2. Therefore, we decided to estimate the cooperativity as that of the dimeric tetR transcriptional regulator. TABLE Parameter Value Source INSERT PHENOMENOLOGICAL TF The reaction between Pyocyanin and soxR was modelled as a fast, reversible reaction with the following equation, where K represents the reciprocal of the concentration at half maximal GFP, and m is a measure of the gradient of the resulting sigmoid. We initially computed the ODE model using both Euler’s method and Matlab’s built-in suite of ODE solvers in order to check that the automated solvers were returning valid results. Once this had been confirmed, we switched to exclusively using the ODE solvers in order to increase model efficiency and perform feature extraction (such as maximum GFP expression and points of inflection). INSERT DIAGRAM COMPARING EULER TO ODE
We assumed that the oxidation of soxR by Pyocyanin was fast compared to the timescales of the biological processes in the model, and therefore modelled the quantities of oxidised and reduced soxR as instantaneously created from the total available pool of soxR, splitting the pool in the following way. INSERT EQUATION FOR PYOCYANIN PROPORTIONS We also assumed that the cells had been left to grow to steady state before the simulation, and calculated the steady state concentrations of soxR and soxR mRNA to use as the initial conditions for the respective ODEs.
Comsol Multiphysics
In order to better understand the system we are engineering, we constructed a deterministic transcription-translation model of the gene circuit using ordinary differential equations in Matlab. This model allowed us to explore the dynamics of the system, to better understand the parameters that the model would be sensitive to, to use curve fitting tools to estimate unknown parameters, and to check whether the constructs were behaving as expected.
Using parameters from literature, we constructed a model consisting of 4 coupled ODEs: one for transcription of soxR, one for translation of soxR, one for transcription of GFP, and one for translation of GFP. SoxR production was modelled as constitutive, as there was insufficient evidence in literature to support feedback control of the transcription of soxR. The binding of soxR to its promoter was modelled using the phenomenological transfer function given below. There was no clear data to be found regarding the cooperativity of the soxR transcription factor; however, it is known that soxR forms dimers, so the maximum cooperativity is constrained at 2. Therefore, we decided to estimate the cooperativity as that of the dimeric tetR transcriptional regulator. TABLE Parameter Value Source INSERT PHENOMENOLOGICAL TF The reaction between Pyocyanin and soxR was modelled as a fast, reversible reaction with the following equation, where K represents the reciprocal of the concentration at half maximal GFP, and m is a measure of the gradient of the resulting sigmoid. We initially computed the ODE model using both Euler’s method and Matlab’s built-in suite of ODE solvers in order to check that the automated solvers were returning valid results. Once this had been confirmed, we switched to exclusively using the ODE solvers in order to increase model efficiency and perform feature extraction (such as maximum GFP expression and points of inflection). INSERT DIAGRAM COMPARING EULER TO ODE
We assumed that the oxidation of soxR by Pyocyanin was fast compared to the timescales of the biological processes in the model, and therefore modelled the quantities of oxidised and reduced soxR as instantaneously created from the total available pool of soxR, splitting the pool in the following way. INSERT EQUATION FOR PYOCYANIN PROPORTIONS We also assumed that the cells had been left to grow to steady state before the simulation, and calculated the steady state concentrations of soxR and soxR mRNA to use as the initial conditions for the respective ODEs.
Comsol Multiphysics
For every model we wished to fit to data, we developed an objective function G of the model F, and the data Y. We then used optimisation algorithms to minimise this function G, by altering values of the parameters params. Insert: min(G = (Y - F(data,params))^2) Caption: Curve fitting regime
UMost notably, this regime was used to estimate parameters for the reaction between Pyocyanin and SoxR, which were unknown, as well as the basal expression, the GFP:RFU scale of the measurements, and the degradation rate of GFP, given the transcription/translation rates sourced from literature. In order to utilise the entire available dataset, rather than merely the steady state GFP values for each Pyocyanin concentration, we constructed an objective function that could generate ODE plots. The reasoning behind using this more complicated strategy is as follows. The steady state GFP reading can be computed using equation [INSERT EQUATION NO] below. It can be seen that there is orthogonality in all of the parameters in table INSERT TABLE NO, besides the Scale and the degradation rate of GFP; a twofold increase in the Scale parameter can be compensated for by a twofold increase in the GFP degradation rate. However, the same is not true in the time-dependent model. It can be seen from figure INSERT FIG NO that a twofold increase or decrease in these parameters leads to a differing transient response, and the same steady state response. Thus, an objective function that utilises transient data can fit both the GFP degradation rate and the GFP:RFU scale simultaneously. INSERT equation: steady state GFP INSERT table: fitted params INSERT fig: Orthogonality
Convergence of optimisation algorithms is of vital importance, and plotting the norm of the residuals (ie the dataset-wide error of the model) can be qualitatively very informative. An optimisation problem can be said to have converged if the norm of the residuals tends towards a certain value. Thus, it can be seen from figure [Convergenceofc1amaster] that convergence is obtained, and from figure [fit25] that a reasonable fit of the data is achieved.
Analytical Modelling Guide
