Team:HKUST/Model

ShellBi

Project
Description Design Engineering Proof of Concept Implementation Safety Contribution Notebook
Wet Lab
Experiments Results
Dry Lab
Model
Reach
Human Practice Collaboration Communication Education Sustainable
Parts
Basic Parts Composite Parts
Team
Team Attributions

Model

Assumptions and Simulation

Overview

The core mechanism of our biosensor involves two main parts –– the antagonism between Paralytic Shellfish Toxins (PST) and veratridine (VTD) on ion influx, and the Two-Component System (TCS) in E. coli. The antagonism of PST and VTD on channel proteins detects the presence of the toxins and induces an ion influx into the cell. TCS is used to express green or red fluorescent proteins (GFP/RFP) according to the osmolarity change.

To improve the performance of the TCS, we have incorporated antisense RNA as well as an engineered scaffold protein with a feedback loop. The antisense works at the frontend to reduce leaky false expression of GFP or RFP, while the scaffold makes the backend TCS response more sensitive and rapid.

Figure 1: circuit diagram of our engineered TCS

To search for evidence for the antagonism, we have tried molecular docking of PST and VTD on an ion channel. Due to the complexity of the circuit, we first investigated the behaviour of each individual part separately, and then we integrated each module to simulate the entire circuit. We adopted a deterministic approach using ordinary differential equations to model the kinetics for each biochemical interaction. The constants in our model are fitted with experimental results from the wet lab and literature. For details of the mathematical equations, please refer to the supplementary information.

Supplementary Appendix

In addition, we have done sensitivity analysis to identify the important parameters in our system. Our goal is to find ways to optimise the performance of our biosensor by fine-tuning each part in the system, especially the antisense and the scaffold.

In summary, our goal for modeling is:

• ❒ Predict potential target for the antagonism
• ❒ Describe the behaviour of the engineered TCS system
• ❒ Highlight the importance of each extra component in the engineered TCS system
• ❒ Optimise the system with exchangeable parts
• ❒ Use the optimised model to predict the situation when saxitoxin is present

Assumptions of Our Model

• 1. We consider EnvZ-OmpR Two Component System is alone and the intermediates do not have any reactions with other systems
• 2. OmpRP only activates pOmpC and pOmpF
• 3. All elementary reactions follows mass-action kinetics
• 4. All osmolytes have uniform effect on altering autophosphorylation rate
• 5. The concentration of EnvZ should be lower than concentration of OmpR
• 6. The tolerance of salt by E. coli is negligible
• 7. Adenosine triphosphate (ATP), which is the phosphate donor, is always in excess

PST effect on Ion Channel

Molecular docking is performed using Autodock vina [1]and CB-dock [2]. The data for the ligand saxitoxin and veratridine are from pubchem while the data for Nav1.7, MscS(open) and Mscs(close) are from RCSB-PDB database. The protein ids are 6j8g, 6vyl and 6vyk respectively.

After random testing, we found that E. coli small-conductance mechanosensitive channels (MscS) have a similar binding score with human voltage-gated sodium channel(VGSC) Nav1.7 when docked with saxitoxin (STX) and veratridine (VTD). For STX binding onto Nav1.7, the vina scores are around -7.6 and 6.4 in autodock vina. The results are demonstrated in Fig(2) and Fig(3)

Figure 2: Autodock vina results for STX and MscS

Figure 3: Autodock vina results for STX and Nav1.7

From the results, the magnitude of the score for STX and MscS is actually higher than that for STX and Nav1.7. Since Nav1.7 is a characterised protein that is known to be inhibited by STX, the higher magnitude score suggests that there are interactions between STX and MscS. To verify it, CB-dock, another docking software, is used to double verify the results and provide easy visualization of the ligand position on the receptor. In addition, MscS interaction with VTD is also tested by the two softwares.

For the results of Autodock vina with veratridine and MscS and Nav1.7:

Figure 4: Autodock vina results for VTD and MscS

Figure 5: Autodock vina results for VTD and Nav1.7

From the results of autodock vina, it demonstrates two important conclusions about the channel. First, VTD binds to MscS more tightly than STX. This is explainable, because in VGSCs, the binding of STX is reversible but the binding of VTD is not. Secondly, the score for binding of VTD to MscS and Nav1.7 is quite similar in the range of 9, suggesting that VTD can also bind well to MscS.

In molecular docking, the score is not the only factor to show that the binding is reasonable. We need to study the result of the docking as well to verify the results. In light of that, we used the same input files with autodock vina to an online docking server CB-dock as it has included molecular visualization for different docking results.

For the CB-dock for STX and MscS and Nav1.7

Figure 6: CB-dock results for STX and Nav1.7 right: STX + MscS

For molecular visualization:

Figure 7: CB-dock visualization of STX binding into Nav1.7

Figure 8: CB-dock visualization of STX binding into MscS

The binding scores comparison is quite similar to autodock vina, so the same conclusions apply.

From the visualization, STX did bind to the desired location in Nav1.7, which is around the p-loop and the transmembrane region of the pore region. Therefore, the score is acceptable as a reference to determine the binding affinity of STX to other proteins. For Mscs, STX binds in the pore region in the channel, which suggests that STX once binds to Mscs, it could potentially inhibit the channel ion influx, which occurs when STX binds to Nav1.7. Therefore, STX has a high potential to inhibit Mscs ion influx during osmoregulation.

For binding veratridine and Nav1.7 and MscS and Nav1.7 in CB-dock:

Figure 9: CB-dock results for VTD and Nav1.7

Figure 10: CB-dock results for VTD and MscS

Figure 11: CB-dock visualization for VTD and Nav1.7

Figure 12: CB-dock visualization for VTD and MscS

The binding scores comparison are also quite similar to autodock vina, so the same conclusions apply.

From the visualization, VTD binds to the transmembrane region of Nav1.7 and therefore matches the ideal result of the VTD binding site in Nav1.7. For the binding of VTD to MscS, the binding affinity is strong, as denoted by the docking score. In addition, the docking position of VTD is in the transmembrane region of the protein. Since this is a mechanosensitive region, extra binding on the transmembrane region potentially induces gating changes, which also occurs in Nav1.7.

To solidify our findings, extra docking using CB dock on the close confirmation of Mscs are docked with STX and VTD.

The binding score for STX & VTD to close Mscs:

Figure 13: CB-dock results for STX and closed MscS

Figure 14: CB-dock results for VTD and closed MscS

From the binding scores, since the channel is now closed, binding affinity for both STX and VTD should drop. The binding score matches expectations. Besides, the score of VTD is still higher than the score of STX, so the previous conclusion also stands.

Visualizations of the bindings:

Figure 15: CB-dock visualization for STX and closed Mscs

Figure 16: CB-dock visualization for VTD and closed Mscs

From the binding sites, STX again binds onto the pore region of Mscs, indicating its potential function as an inhibitor to the channel. For VTD, it also binds to the transmembrane region of the channel, suggesting that VTD could potentially bind to MscS when the channel is close, supporting the previous hypothesis of VTD changing the gating of the channel.

As a result, we believed that MscS is a potential target for the antagonism. Combined with the fact that Mscs is one of the most abundant channels for osmoregulation in E. coli [3], we believe our system will definitely work.

Two Component System

The two components in the TCS in our case are EnvZ and OmpR. EnvZ senses intracellular osmolarity, which causes downstream phosphorylation of OmpR [4]. EnvZ senses internal but not external osmolarity [5]. Phosphorylated OmpR can act as a transcription factor for the promoters PompC and PompF. Under PompC and PompF we express GFP and RFP respectively.

The osmolarity alters the rate of EnvZ’s autophosphorylation (kap), which refers to the rate of EnvZ gaining a phosphate group [6]. However, from literature about PST and VTD antagonism, only intracellular sodium level was measured [7]. Therefore we have used [Na+] instead of osmolarity as our input. By fitting to experimental data, the effect of [Na+] on kap can be estimated in the form of a second-order Hill’s equation:

Figure 17: hill's equation for EnvZ autophosphorylation

After obtaining a phosphate group, EnvZ would act as a kinase and transfer it to OmpR. Note that non-phosphorylated EnvZ would act as phosphatase and remove the phosphate group from OmpR~P. For a summary, the model can be represented by the pathway [8]:

Figure 18: Simplified pathway for EnvZ/OmpR system

We have adopted the binding, phosphotransfer, and dissociate interaction as our model. All interactions are considered as first order [8].

Figure 19: pathway for phosphate transfer

The details of the biochemical equations are as follows:

In (5), EnvZ cannot undergo phosphorylation when it is bound to non phosphorylated OmpR, forming a “dead-end complex” [8]. Dissociation is required to allow autophosphorylation of EnvZ.

Phosphorylated OmpR would then dimerise and then bind to DNA to regulate gene expression. There are various models proposed for OmpR dimerisation, but it has been elucidated that the most probable sequence of OmpR mode of action is phosphorylation, dimerisation, and DNA binding [9].

PompC and PompF

The concentration of phosphorylated OmpR dimer (OmpR~P) determines which promoter to be activated. In low osmolarity, a low level of OmpR~P is present, and PompF is preferentially activated. Conversely, in high osmolarity, high levels of OmpR~P correspond to preferential activation of PompC [10].

The mechanism of the reciprocal expression of PompC and PompF involves hierarchical and cooperative binding of OmpR~P. In both PompC and PompF, there are several binding sites for OmpR~P, namely C1, C2, C3, F1, F2, F3, and most importantly, the F4 site. All of these sites, except F4, when bound by OmpR~P, would activate the promoter. Furthermore, as mentioned above, OmpR~P is a dimer, and here each site has two subsites for each monomer of OmpR, termed C1a, C1b, C2a, C2b, and so on, with all b-subsites having higher affinities to OmpR~P than a-subsites except F4 [10].

Figure 20: PompF and PompC binding

The sites have a preferred order to be bound, which is experimentally determined to be F1, C1 > F2, F3 > C2 > C3 [10]. It means that at a low level of OmpR~P, F1 and C1 sites are bound first, and as the OmpR~P level goes up, F2 and F3 would also be bound by OmpR~P, so on and so forth. This explains that at low OmpR~P, they bind to F1-F3, leaving C2 and C3 unbound, and thus PompF is activated but not PompC.

The remaining F4 site, which is not in the above hierarchy, has the opposite effect of deactivating the transcription. At low OmpR~P, F4 is not bound, leaving PompF to be activated. However at high OmpR~P, all F1-F4 are bound, and F4 restricts the transcription from PompF [10]. As a result, at high OmpR~P, the binding of OmpR~P on the F4 site represses PompF, while the binding of OmpR~P to C1-C3 allows PompC to be activated.

Furthermore, the 3 activating OmpR~P cooperatively bind to both C1-3 and F1-3. Binding to the F1 site helps increase the OmpR~P binding affinity of the neighbouring F2, and F2 would increase F3 affinity, and the same applies for C1 to C2 and C2 to C3 [10]. As mentioned above, there are inconsistencies about the order of OmpR phosphorylation, dimerisation and DNA binding between different studies. It is proposed that binding of monomer OmpR to the b-subsite cooperatively helps another monomer to bind to the a-subsite [10]. However, since there is stronger evidence that Omp~P dimerises before binding to the DNA, we do not consider cooperative binding between subsites [9].

Combining the effect of hierarchical and cooperative binding, we propose that the expression level of PompC and PompF could be described by Hill’s equations:

Both Hill’s coefficients nC and nF have an theoretical upper bound of 3, which is the unrealistic case of maximum cooperativity of the C1-3 and F1-3 sites. The F4 site is a negative form of Hill’s function with no cooperativity and it is independent of the other F sites.

The model we found in literature uses first-order Hill’s function only [11], but we found that non-trivial Hill’s coefficient matches experimental data better, by having a sigmoidal curve for PompC [12]. By fitting the wet lab data, we found that it estimates the behaviour well when both Hill’s coefficients are 2.

Finally, for engineering, the promoters could be tuned by mutating the ab-subsites. Recall that the b-subsites bind more tightly to OmpR~P than a-subsites. It is found that by mutating the F1 site into different combinations of F1a and F1b, the binding affinity to OmpR~P was observed to be F1aF1a < F1aF1b (wildtype) < F1bF1b, and the same applies to other sites [10]. The binding affinity could be modelled by altering the dissociation constants KC and KF. Hill’s coefficients nC and nF are approximated to be constants, as the cooperative binding of OmpR is due to the interactions of protein domains and it should be independent of the promoter DNA sequence.

Antisense

After consulting Professor Aifeng Li from the Ocean University of China, and also from SDU Denmark and University of Indonesia iGEM teams, we learnt that there are leaky expressions from PompC and PompF. Thus antisense RNA is used to repress the leaky expression of GFP and RFP when they should not be expressed.

An antisense RNA is complementary to a fragment of mRNA. It works by complementary pairing to a coding mRNA to block translation. In our circuit, we put the antisense of GFP (asGFP) under the promoter PompF, as opposed to its target GFP, and the antisense of RFP (asRFP) under PompC. Please refer to design for more details.

Figure 21: Antisense pathway

We modelled the production of GFP and RFP in a two-step process of transcription and then translation. The mechanism of antisense is modelled by considering the binding and degradation. First is the formation of mRNA-antisense complexes, stated in (8) and (11), with the degradation of both mRNA and antisense in their free unbound form (not shown in the equations below), and lastly, the degradation of both mRNA and antisense individually in their complex form, shown in (9), (10), (12) and (13) [13].

Furthermore, we have included the GFP and RFP maturation by protein folding, which can capture the time delay for the colour to be observable [13].

Scaffold

We had an interview with a local oyster farmer to learn more about the feasibility of our product as well as the needs of our consumers. The oyster farmers expect the biosensor to be fast and sensitive. Therefore, we have modified the TCS by introducing a third protein as a scaffold (Sc) to boost the sensitivity and rapidity of the system.

The EnvZ is modified by adding an SH3 domain (SH3dom) as well as inserting an SH3 ligand (SH3pep), such that without the scaffold, the SH3 domain would bind to the ligand and block phosphorylation of OmpR, which is called autoinhibition [14].

Figure 22: autoinhibition and binding of scaffold to EnvZ

For the OmpR, we added a monomer of a leucine zipper (LZ), which pairs with another monomer of the same leucine zipper. The third player, the scaffold, is a chimera of SH3 ligand and the complement monomer of the leucine zipper, so that it can bridge from EnvZ to OmpR and unlock EnvZ from its autoinhibition.

Figure 23: Effect of scaffold to EnvZ/OmpR system

Based on the model of the TCS, we have modified the kinetics by including the formation and dissociation of all possible combinations of dimers involving the scaffold by (19)-(22) and trimers by (23)-(28) [15]. EnvZ autoinhibition is written as (17). We assume binding and dissociating of SH3 and LZ are invariant despite the species involved in the reactions.

According to the literature, both free OmpR~P and scaffold bound OmpR~P are capable of being transcription factors [15]. In our case, since EnvZ is localised on the membrane, the EnvZ-Scaffold-OmpR~P complex is not considered to be a transcription factor. This is expressed in (37), and this would be the input for PompC and PompF in (7) and (8).

From the literature, there are 3 candidates for the SH3 ligands and 2 pairs of leucine zippers with different dissociation constants [14]. This will be discussed in the Optimisation section.

Optimisation section

A feedback loop is formed by putting the scaffold under PompC, which will be presented in the next section. For more details on the feedback mechanism, please read the design.

Design page

Model Calibration and Simulation

Therefore, with all the chemical reactions, we construct a mathematical model based on the following equations:

Detailed Equations

More in Appendix

Fitting wet lab data

From our wet lab, we obtained data for intracellular osmolarity versus extracellular osmolarity and tested Two Component Systems with eGFP and mRFP1 as reporters in E. coli BL21. The results are as follows:

Figure 24: Wet lab data for eGFP intensity at different osmolarity

Figure 25: Wet lab data for mRFP1 intensity at different osmolarity

Detailed results can be found below

Results page

Combined with the known constants from literature, the list of rate constants and their corresponding values are shown below

Supplementary Appendix

The fitting results are as follows

Figure 26: Fitting results with our TCS basal model to wet lab data

By preliminary analysis, both curves are nicely fitted with wet lab data, so the trend of our model matches in vivo experiments. However, not the whole range of intracellular osmolarity can be captured from the above figure. In addition, The promoter strength of GFP is much stronger than RFPs. Therefore, in the following discussion the promoter strength constant is tuned down for analysis.

Since we have the parameters and a verified model now, we can predict the behaviour of wild type TCS, wild type TCS with antisense, engineered TCS without antisense and finally engineered TCS with both antisense and scaffold, and compare them to highlight the importance of antisense and scaffold in a graphical format.

Basal Wild Type TCS

For wild type TCS, the simulation graph of RFP from PompF and GFP from OmpC against intracellular osmolarity is shown below. It can be seen that RFP from the promoter PompF forms a spike curve at around 18 mOsm/kg for its highest amplitude while GFP from the promoter PompC forms a sigmoidal curve and levels after 50 mOsm/kg. This shape can be explained by below.

Hierarchical and cooperative binding of OmpR~P

When intracellular osmolarity is low, PompF is activated while there are only minimal expressions for PompC. However, at higher osmolarity, F4 site on PompF is bounded and inhibits the expression of PompF, so the level of expression of RFP drops. At the same time, PompC is activated at high osmolarity and more GFP is expressed. At higher osmolarity, since all OmpR has been phosphorylated, the curve levels off. However, the system itself is not bistable as at 20 mOsm/kg, GFP expression activates but RFP is still expressing, so two colours of fluorescent proteins are expressed and harder to observe. Besides, the sensitivity of the system, which is denoted by the intersection point, is not tunable. This makes the system unsuitable since the results are not obvious and insensitive.

However, the system itself is not bistable as at 20 mOsm/kg, GFP expression activates but RFP is still expressing, so two colours of fluorescent proteins overlap and become harder to be analysed. Besides, the sensitivity of the system, which is denoted by the intersection point, is not tunable. This makes the system unsuitable for a biosensor since the results are ambiguous.

Figure 27: TCS basal simulation results

Scaffold

With the addition of scaffold to TCS basal, PompC activates minimally at lower osmolarity. When the osmolarity rises, it triggers a positive feedback loop such that when scaffold unlocks EnvZ from autoinhibition, the concentration of OmpR~P rises rapidly which triggers more scaffold to be expressed until it reaches the point where the rate of scaffold production equals to its rate of degradation. As a result, the presence of the scaffold reduces RFP signal from PompF when the expression of PompC is greater than PompF, reducing the overlapping of positive and negative signal, which can be confusing and may lead to huge monetary losses to farmers, referring to the interview of the oyster farmer.

Interview of the oyster farmer

Besides, the positive feedback loop created by expressing scaffold on PompC allows even smaller expression of GFP from PompC to amplify and reach its maximum, creating a step-function curve when PompC is stronger than PompF. In addition, the bandwidth of the RFP signal is dependent on the strength of expression of the scaffold, which can be easily tuned. More about this can be found in the optimisation section.

Optimisation section

Figure 28: TCS with scaffold simulation result

In order to compare with TCS basal, we overlay the two graphs into one for analysis, which is shown below. For their amplitude, which is GFP and RFP expression level, both graphs reach identical peaks. This indicates that the scaffold does not affect the expression of the reporters, which suits our interest as we don’t want to weaken our reporters. For their intersection point, to make the intersection more clear, the intersection point for TCS basal is coloured in magenta while the intersection point for TCS with scaffold is coloured in cyan. From the simulation graph below, when the system is equipped with the scaffold, the intersection point shifts left due to the positive feedback loop. However, the concentration of fluorescent protein remains unchanged at the two intersection points. This suggests that the shift in the scaffold is linear, thus making it easier to engineer. At the same time, the curve bandwidth of RFP is much sharper than TCS basal, so once the the result reaches the threshold, green signal is strongly expressed so false positive rate can be dramatically reduced.

Figure 29: comparison between TCS scaffold and TCS basal

Antisense

In our design, we added complementary antisense to PompF and PompC, so depending on the expression level of PompF and PompC, the one with a higher expression should suppress the other fluorescent protein, thus creating a bistable-like response depending on whether the osmolarity is higher than the intersection point or lower than the intersection point. As shown in the simulation graph below, although the curve shape remains unchanged, the expression level of both GFP and RFP is lowered, as denoted by the y-axis. The expression level can be tuned as well.

Optimisation section

Figure 30: TCS with antisense simulation result

In order to compare with TCS basal, we overlay the two graphs into one for analysis, which is shown below. From the simulation graph, the intersection point for TCS basal is coloured in magenta and the intersection point for TCS with antisense is coloured in black. By comparison, the amplitude of both GFP and RFP decreases a lot, but the region where there are two signals together are reduced as well due to antisense suppressing the weaker fluorescent protein into negligible amounts. In addition, the intersection point does not shift laterally but only vertically, so it can be easily engineered as well.

Figure 31: comparison between TCS with antisense and TCS basal

TCS with scaffold & antisense

Combining the effect of antisense and scaffold, the simulation is shown in the graph below. Both effects of antisense and scaffold is incorporated into the final circuit, where the scaffold can create a positive feedback loop to produce a sharp step-function curve when the osmolarity is high, while the antisense suppresses complementary signal, so the low expression of PompC and relatively higher PompF expression at low osmolarity suppress each other. However, at high osmolarity, PompF is already inhibited by itself while PompC expresses at a high level, so the effect of antisense is not significant at high osmolarity.

Figure 32: TCS with scaffold and antisense simulation result

To compare it with the TCS basal, the simulation graphs are overlaid again. As mentioned in section TCS with antisense and TCS with scaffold, the effects of antisense can be simplified into compression of the curve in y-direction, while the effects of the scaffold can be simplified into compression towards the x-axis. As denoted by the graph below, the blue intersection point of TCS with scaffold and antisense shifted both downwards and leftwards relative to the magenta TCS basal intersection point, thus proving their effects respectively.

Figure 33: comparison between TCS basal and TCS with scaffold and antisense

In summary, the graph of our engineered TCS is shown below. The result of our engineered TCS reduces a mix of red and green coloursto be observed while the false-positive rate and sensitivity of the system can be controlled at the same time.

Figure 34: Final simulation result for our engineered TCS system

Scaffold only vs Scaffold and Antisense Component-wise, for the function of antisense, when comparing TCS with scaffold and our engineered TCS, the effects of antisense RNA, which is suppressing the complementary signal, are demonstrated in the graph below. The magenta intersection point for TCS with scaffold lies on the same x-axis as the black intersection point of our engineered TCS while it only lowers its level. Therefore, as mentioned above, the effect of antisense is observable in our engineered TCS.

Figure 35: Comparison between our engineered TCS system and TCS with scaffold

Antisense only vs Scaffold and Antisense For the function of the scaffold, when comparing TCS with antisense and our engineered TCS, the effects of the scaffold, which is reducing the bandwidth of the PompF curve, is shown in the graph below. In terms of the intersection points, the black intersection point of our engineered TCS shifts left and down relative to the magenta intersection point of TCS with antisense. The horizontal shift is explained as the scaffold induces PompC expression. However, this induction also triggers antisense production, so RFP antisense expression at low osmolarity becomes higher, which explains the vertical depression of the intersection point. In addition, the PompC curve exhibits a digital signal, in which there is a steady intensity when the osmolarity exceeds the intersection. This allows the GFP signal to be amplified once the intersection point is reached, making it user-friendly and the signal can be distinguished quickly compared to an analog curve.

Figure 36: Comparison between our engineered TCS system and TCS with antisense

Therefore, our engineered TCS is a user-friendly and sensitive biosensor that allows accurate detection of intracellular osmolarity changes.

Sensitivity Analysis

In the complete TCS, ordinary differential equations (ODE) are utilized to model the pathway as a mathematical expression of the system. There are several parameters in the system that cannot be accurately determined and the parameters in other experiments cannot be directly used in our system since they are based on different experimental conditions. Therefore, sensitivity analysis (SA) needs to be carried out in order to identify the influential parameters and to help us to optimize the model.

Sensitivity analysis is the study of how the uncertainty in the output of a mathematical model or system (digital or other) is divided and allocated to different sources of uncertainty in its input [16, 17]. Sensitivity analysis is mainly divided into two categories: local SA and global SA. Global sensitivity analysis is the process of assigning the uncertainty of the output to the uncertainty of each input factor in its entire range of interest, while local sensitivity analysis evaluates the local impact of changes in input factors on the model response by focusing on the sensitivity around a set of factor values [18, 19]. In our system, we use global SA since all input parameters will change at the same time and sensitivity will be evaluated across the entire range of each input. Therefore, Sobol analysis was adopted for global analysis. While Morris’s method is less complex but doesn’t provide much information about specific interactions between variables, Sobol analysis requires demanding computations but it is more accurate [20], we decide to use Sobol analysis to initially screen out the important parameters and then modify them to optimize the system.

Sobol analysis, also known as variance-based sensitivity analysis, provides first-order and total indices, which respectively describe the fraction of output variance contributed by each factor on its own, and by the sum of first-order and all higher-order interactions for each factor.[20, 21]

For the Sobol method, we treat the model as a function

where X is a vector of d uncertain model inputs as below.

Decompose f(X) in the following way,

The total unconditional variance can be defined as

where Ω^d represents the d-dimensional unit hyperspace. In order to get the partial variances which are the components of the total variance decomposition, we need to compute each of the terms in equation (39) as

Then we can determine the variance contributions to the total output variance and interactions between parameters. The contributions are written as Sobol sensitivity indices as below,

Since there are several parameters for the reactions we need to determine in the TCS and the rough range of them were obtained from experimental data of previous research, we first applied Sobol sensitivity analysis to all the parameters in order to find out the parameters that have a great contribution to the output. As for the Sobol analysis, we studied the method and code from the Thessaloniki team in 2018, which used SALib Library for sensitivity analysis.

After computation, the figure below shows the first order (S1) and total order (ST) sensitivity indices for all the parameters. Orange bars represent the first-order indices and the total order sensitivity indices which are larger than one are shown in steel-blue and we treat them as important total order sensitivity indices (iST), while ST less than one are shown in light sky blue. As we can see from the figure, the first order indices are quite small compared to the total order indices, which means that the interactions between the parameters play an important role in the system instead of separate parameters, and the system has a high complexity since there are eight parameters whose ST are larger than one. This result confirms that using global sensitivity analysis instead of local analysis is more suitable for our complex system. Therefore we successfully picked out 8 important parameters (kd2, kb2, kph, kd3, kb3, kbLZ, kdLZ and kbSH3) among 23 parameters so that the difficulty of the evaluation can be reduced.

Figure 37: Sobol analysis results for all parameters

In order to further confirm the importance and influence of the selected parameters, total order sensitivity indices are computed under different sample sizes. The accuracy of the sensitivity indices increases when the sample size increases. As is shown in the figure below, the parameters fluctuate greatly in the early stage, as the number of samples increases, the values tend to be stable and have no significant difference from each other.

Figure 38: Sensitivity indices under different sample sizes

Using Sobol analysis, the 3 most important parameters are kph (rate constant for the dephosphorylation of OmpR by EnvZ), kbLZ (rate constant for the binding of leucine zipper) and kbSH3 (rate constant for the binding of SH3). However, kph is inherent to EnvZ so only kbLZ and kbSH3 will be varied to optimize the system.

Optimisation

For our engineered system, some parameters are fitted and therefore may not be optimised and suits our interest to provide a user-friendly and versatile biosensor. Therefore, combined with the results from sensitivity analysis, optimization is done on the exchangeable parts.

There are primarily 3 components that are optimized, where are the scaffold, antisense, and the PompC and PompF promoters. Note that the amplitude is significantly different from the original simulation result as some of the parameters are tuned down to enhence the visualisation of the effect of tuning the parts.

Scaffold optimisation

As the scaffold is very flexible, the scaffold binding affinity, degradation and expression level can all be tuned.

From sensitivity analysis, kbSH3 and kdLZ are important parameters mathematically. In a physical sense, kbSH3 is the binding rate of the scaffold to EnvZ while kdLZ is the binding rate of the scaffold to OmpR. From the literature, there are different parts with different binding affinities [14]. Since we now have a model to predict its effects, they are numerically varied to explore the effects of different binding affinities and choosing the preferential part.

For kbSH3

To test for kbSH3 only, the binding affinity of the scaffold to OmpR is set very high to allow kbSH3 to be a rate-determining step while the autoinhibition of EnvZ is also set very high to amplify the effect.

The results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. For GFP, the effect of kbSH3 is not very significant since the scaffold is already in excess. However, for RFP, since the scaffold expression is low, the scaffold needs to overcome autoinhibition of EnvZ so the binding rate kbSH3 must be higher than the autoinhibition binding rate. From the right plot, the graph demonstrates the effects of autoinhibition. The presence of autoinhibition contributes to the sudden drop in RFP concentration at 10 mOsm/kg when kbSH3 is lower than 10, which is the value of EnvZ autoinhibition binding affinity.

Therefore, to overcome the autoinhibition, the binding affinity of the scaffold must be designed to be higher than EnvZ autoinhibition.

Figure 39: GFP concentration for variations in kbSH3

Figure 40: RFP concentration for variations in kbSH3

For kbLZ

To test for kbLZ only, the binding affinity of the scaffold to EnvZ is set very high to allow kbLZ to be a rate-determining step.

The results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. For better visualisation, the graph's x-axis representing osmolarity is shown in the opposite direction. In general, when kbLZ is higher, more OmpR will bind to the scaffold. However, this is not necessarily proportional to the scaffold horizontal shift because when the scaffold binding affinity to OmpR is smaller, the scaffold bound to OmpR~P can be reused each time to form new scaffold and OmpR complexes. We found that the recycling effect is more significant than the abundance of the scaffold, so from both RFP and GFP, low kbLZ actually shifts the curve leftwards and reduces the bandwidth.

Therefore, to tune the sensitivity of the system, kbLZ can be tuned so that our TCS system can detect PSTs and other neurotoxins at a broader range of induced osmolarity. In terms of parts, we should use the LZ with a lower binding affinity, which increases kbLZ to increase the bandwidth of the positive region.

Figure 41: GFP concentration for variations in kbLZ

Figure 42: RFP concentration for variations in kbLZ

For scaffold degradation rate

To test for the significance of degradation rate, the expression rate constant for the scaffold is in the middle of the tested range of degradation rate constant to amplify the effect.

The results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. For better visualisation, the graph's x-axis representing osmolarity is shown in the opposite direction. When the degradation rate is higher, the steady-state level of the scaffold decreases, so the effects of scaffolding to amplify PompC decrease. As a result, the curves shift right when the degradation rate is increased.

Compared to kbLZ, although the effect is similar, the effect by degradation is more straightforward so it is more advisable to tune the degradation rather than kbLZ. Note that the tuning effect is saturated once it reaches a certain level, so it is neglected in the sensitivity analysis despite having an effect here.

Figure 43: GFP concentration for variations in scaffold degradation

Figure 44: RFP concentration for variations in scaffold degradation

For scaffold expression rate

To test for the significance of expression rate, the degradation rate constant for the scaffold is in the middle of the tested range of expression rate constant to amplify the effect.

The results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. For better visualisation, the graph's x-axis representing osmolarity is shown in the opposite direction. When the expression rate constant of the scaffold increases, the curve shifts left because higher scaffold production can amplify the PompC and therefore GFP expression will be higher, so the curve shifts to the left in the case of high expression rate constant.

Therefore, since it exhibits reciprocal effect with degradation rate, they can be used to tune the system for sensitivity. Note that the tuning effect is also saturated once it reaches a certain level, so it is ignored by the sensitivity analysis despite showing an effect here.

Figure 45: GFP concentration for variations in scaffold expression

Figure 46: RFP concentration for variations in scaffold expression

Antisense optimisation

Antisense specific to GFP and RFP can be designed such that the expression, degradation and binding affinity can be varied based on the length or base pair modification of antisense fragments. Therefore, we try to vary their values to observe and optimise the antisense parts.

For binding affinity of GFP antisense & RFP antisense

To test for the binding affinity of RFP antisense and GFP antisense to their targets, the degradation of RFP antisense and GFP antisense is set at a minimal number while expression of both antisense increases so the effect of varying binding affinity can be amplified.

For binding rate of RFP antisense, the results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. GFP is also considered here to see whether tuning RFP antisense also has an effect on GFP. From the graphs below, GFP concentration exhibits minimal effect. However, for RFP, although the binding affinity of RFP antisense is increased, the region where the RFP antisense is not expected to be expressed, in this case as high osmolarity, already has a low amount of RFP being expressed due to inhibition of PompF. The feedback loop formed by the scaffold further provides a synergistic effect. Since the binding affinity only allows the binded mRNA RFP to not be released back into its free form, the strength of antisense is masked by the other two mechanisms, so its effect is insignificant.

Figure 47: GFP concentration for variations in binding affinity of RFP antisense

Figure 48: RFP concentration for variations in binding affinity of RFP antisense

For binding rate for GFP antisense, the results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. RFP is also considered here to observe the potential effect of GFP antisense. Similarly, the effects of antisense binding affinity is not significant because once the opposite promoter is expressed, the complementary reporter is automatically preferentially expressed and PompC at low osmolarity is not significant and binding affinity itself only locks a minimal amount of mRNA GFP, so the effect of antisense is not significant.

Figure 49: GFP concentration for variations in binding affinity of GFP antisense

Figure 50: RFP concentration for variations in binding affinity of GFP antisense

For the expression of GFP antisense & RFP antisense

To test for the expression of GFP antisense and RFP antisense, the degradation of GFP antisense and RFP antisense is set at a minimal number so the effect of varying binding affinity can be amplified.

For the expression of GFP antisense, the results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. RFP is also considered here to observe the potential effect of GFP antisense. From the graphs below, RFP concentration remains the same but the slope of GFP concentration decreases once it reaches the intersection point. Unlike the binding affinity, the surge of GFP antisense can bind to GFP and suppress some GFP signals even when PompF is relatively weaker. Thus there is some suppression in low osmolarity but no observable suppression at higher osmolarity because PompF is further weakened.

In ideal situation, it is desired to have GFP signal immediately rise up to its maximum after passing through the intersection point between PompF and PompC, but at the same time the signal should not overlap much when the osmolarity is near the intersection point, so a trade-off between those two should be evaluated.

Figure 51: GFP concentration for variations in expression of GFP antisense

Figure 52: RFP concentration for variations in expression of GFP antisense

For the expression of RFP antisense, the results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. GFP is also considered here to observe the potential effect of RFP antisense. From the graphs above, GFP concentration remains the same but RFP concentration varies. As there are minimal but significant expressions of PompC at low osmolarity, RFP antisense is expressed and the surge can actually lead to significant suppression of RFP at low osmolarity.

For optimisation, there is a trade-off between the RFP signal bandwidth and intensity, in which the presence of RFP antisense can filter out the lower expression range, but the RFP signal intensity drops.

Figure 53: GFP concentration for variations in expression of RFP antisense

Figure 54: RFP concentration for variations in expression of RFP antisense

For the degradation of GFP antisense & RFP antisense

To test for the degradation of GFP antisense and RFP antisense, the expression of GFP antisense and RFP antisense is set at a high number so the effect of varying binding affinity can be amplified.

For the degradation of GFP antisense, the results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. RFP is also considered here to observe the potential effect of GFP antisense. From the graphs below, GFP antisense degradation exhibits minimal effect on both GFP and RFP concentration because the bounded complex of GFP antisense and mRNA GFP makes the amount of free GFP antisense small, so less free GFP antisense can be degraded even when we add degradation sequences on GFP antisense, thus inducing negligible effect. Note that the purple wall-like representation in the graph of GFP concentration is due to the 0 degradation rate of GFP antisense, so all GFP is suppressed.

Figure 55: GFP concentration for variations in degradation of GFP antisense

Figure 56: RFP concentration for variations in degradation of GFP antisense

For the degradation of RFP antisense, the results are shown in the 3D surface plot below. The GFP concentration is in the left plot while the RFP concentration is in the right plot. GFP is also considered here to observe the potential effect of RFP antisense. From the graphs below, RFP antisense degradation also exhibits very minimal effect on both GFP and RFP concentration because when the expression level of RFP antisense is high, most RFP antisense will bind to mRNA RFP and therefore will be harder to degrade once in a complex form, so a degradation sequence on RFP antisense may also not induce a huge effect.

Figure 57: GFP concentration for variations in degradation of RFP antisense

Figure 58: RFP concentration for variations in degradation of RFP antisense

PompC & PompF optimisation

PompC and PompF can be engineered in a way that the ab-sites can be replaced and their sequential order can be modified to tune the binding affinity of PompC and PompF. For example, it has been tested that putting 2 F1a sites together to form F1aF1a has a weaker binding affinity than wild type F1aF1b sites [10]. Therefore, the promoter itself is tunable. From the equation stated in the below link,

PompC and PompF section

PompC has 1 binding affinity constant while PompF has 2 binding affinity constant, denoted as PompF and PompF4, where PompF region activates it while PompF4 region inhibits it. For simplicity, we only simulate the case when both the GFP and the scaffold are expressed under the same tuned PompC.

PompC

For the variation in the binding rate of PompC (KC), the results of GFP and RFP concentration are shown in the two graphs below, where the left graph represents GFP concentration while the right graph represents RFP concentration. From the simulations, in general, when KC increases, the intersection point shifts right because the drop in binding affinity of the PompC region induces a rise in the concentration required to activate PompC at a higher osmolarity. The curve for RFP is also shifted right with a widened bandwidth because of the scaffold, where higher KC reduces expression of the scaffold and it requires a higher osmolarity to turn on the positive feedback loop to inhibit PompF.

In terms of optimization, as the increase in KC does not saturate at a certain point, this can be the major engineered part to tune the sensitivity of the biosensor.

Figure 59: GFP concentration for variations in binding affinity of PompC

Figure 60: RFP concentration for variations in binding affinity of PompC

PompF

For variations in the binding rate of the activating region of PompF (KF), the results of GFP and RFP concentration are shown in the two graphs below, where the left graph represents GFP concentration while the right graph represents RFP concentration. From the simulations, in general, when KF increases, the bandwidth of RFP curve decreases because although KF should exhibit similar behaviour as KC, since it is not expressing the scaffold, it does not possess the ability to change concentration spike easily. Instead, since KF is not affecting PompC, PompC still functions well at higher osmolarity, causing antisense and scaffold to suppress RFP, cancelling the shifting to the right. In addition, the inhibitory site F4 for PompF is still present and blocks its transcription, so the curve looks like it is becoming smaller in amplitude and width, but in reality it should also shift right, similar to KC.

Figure 61: GFP concentration for variations in binding affinity of PompF

Figure 62: RFP concentration for variations in binding affinity of PompF

To better visualize the effects of variations of PompF KF, the variation is also done on the original TCS with no extra components interfering with the variations. The 3D surface plots are shown below. The left plot represents GFP concentration while the left plot represents RFP concentration. When KF increases, GFP expression remains the same but RFP expression curve shifts right, demonstrating the above explanation.

Figure 63: GFP concentration for variations in binding affinity of PompF in TCS basa

Figure 64: RFP concentration for variations in binding affinity of PompF in TCS basal

PompF4

For variations in the binding rate of the repressing region of PompF4 (KF4), the results of GFP and RFP concentration are shown in the two graphs below, where the left graph represents GFP concentration while the right graph represents RFP concentration. From the simulations, the curve for RFP increases in bandwidth when the KF4 increases because the inhibition by the F4 site only occurs at a higher osmolarity now. While for GFP when there is higher expression of PompF at higher osmolarity, antisense is expressed and thus the concentration of GFP drops and curves downwards.

Thus tuning KF4 can allow the sensitivity to be adjusted with the scaffold activation being unaffected, so that the positive feedback loop also starts at the same osmolarity to create a vertical increase in GFP to its maximum.

Figure 65: GFP concentration for variations in binding affinity of PompF4

Figure 66: RFP concentration for variations in binding affinity of PompF4

To better visualize the effects of variations of PompF KF4, the variation is also done on the original TCS with no extra components interfering with the variations. The 3D surface plots are shown below. The left plot represents RFP concentration while the left plot represents GFP concentration. Similarly, GFP expression is unaffected while RFP expression increases and the bandwidth of the RFP curve increases as more OmpR~P is required to activate KF4 region when the binding affinity is larger.

Figure 67: GFP concentration for variations in binding affinity of PompF4 in TCS basal

Figure 68: RFP concentration for variations in binding affinity of PompF4 in TCS basal

Potential usage and the engineering flexibility in PompF and PompC

Since all the ab-subsites are switchable, it is not necessary to put an inhibitory site on PompF, and PompC can be designed so that it inhibits itself at a certain osmolarity. In other words, their roles can be switched. The scaffold can be expressed under a PompC that differs from the PompC that drives GFP expression. Even a tuned F4 site can be added to the scaffold as a repressor at higher osmolarity to engineer its expression range. Therefore, potentially the scaffold expression can be controlled at high osmolarity, and the PompF can be activated at a higher osmolarity. By altering the sequential order of OmpR~P binding on the separate sites, both promoters can be engineered to perform diverse tasks. For our PST biosensor, we only need to tune KC such that the intersection spot matches the osmolarity when a limit of PST is present in the shellfish samples. More details on implementation.

Implementation page

Conclusion

In our modeling section, we have demonstrated the PSTs antagonism through in silico molecular docking, described the behaviour of our engineered TCS circuit through mathematical modeling, highlighted the importance of each component through comparisons in different simulations and tried to optimise the system through sensitivity analysis and simulations with varying parameters. When all the parameters are tuned and tested, the ideal result is plotted in Fig(69). Overall, we show that our engineered TCS system is highly flexible and versatile bistable-like biosensor for the detection of saxitoxin and potentially other neurotoxins, which we believe can be successfully implemented in real life.

Figure 69: Ideal engineered TCS circuit

Appendix

Phylogenetic tree

Background

As the STX/VTD antagonism has only been tested on cyanobacteria, there is no direct experiment of the antagonism on E coli [7]. At the same time, the cyanobacteria proven to show antagonism, namely Cylindrospermopsis raciborskii, has a long doubling time that is impractical to be engineered within one year of iGEM. Therefore, we set a backup plan of searching for an engineerable cyanobactera strain that has close relation to C. raciborskii, Vibrio fischeri (another bacteria tested for antagonism) or human, so that it might demostrate the antagonism to suit our needs. After meeting with Professor Brett Neilan, the author of the paper testing C. raciborskii, we are advised to build a phylogenetic tree based on ion channels.

Reason to choose MEGA-X

MEGA version X provides a graphical interface and multiple functions such as protein alignment and building phylogenetic trees with different parameters and algorithms [22], which is ideal because it is user-friendly and multiple parameters can be used with each phylogenetic tree.

Detailed process for the phylogenetic tree

All protein sequences are obtained from the NCBI database.


First iteration:
Since from literature calcium, sodium and potassium channels form the precursor for VGSCs [23], we conducted an initial global analysis of 150 randomly sampled proteins related to keywords “potassium”, ''sodium” and ''calcium” in bacteria and use well-characterised eukaryotic VGSCs including Homo sapien and Rattus norvegicus. In the global phylogeny, we noticed a cluster of VGSCs as shown in Figure A1.

Figure A1: MLOriginal1.pdf

We then expanded on this cluster and conducted more bootstrapping on channels in this cluster along with a negative control as shown in Figure A2.

Figure A2: small_aligned_2.pdf

In the verified cluster, it can be seen that except for the negative control, all of the channels, no matter eukaryotic or prokaryotic, characterised VGSCs while no transporters such as antiporters and symporters are not in the same cluster. We can draw a conclusion that bacterial VGSCs and eukaryotic VGSCs share a close common ancestor, so probably bacterial ion channels can be related to the antagonism. However, no significant conclusion on the best fit channel protein can be drawn, so the process cycles into the next iteration.

Second iteration:
In the next cycle, with the suspicion that other bacterial ion channels can provide evolutionary relationships between VGSCs in different bacterial species, we randomly sample bacterial ion channels to verify the previous conclusion. 38 proteins ranging from potassium channels to mechanosensitive channels are used to construct a phylogeny. In the results in Figure A3, we notice the same cluster that appeared in the previous results.

Figure A3: small_aligned_3_extra_3.pdf

We then expanded on this cluster and conducted more bootstrapping on channels in this cluster along with some other negative control as shown in Figure A4.

Figure A4: small_aligned_3.pdf

In the second iteration of the cluster, almost the same cluster is used to perform phylogeny. This indicates that other ion channels share a distant evolutionary path than VGSCs, which indicates that VGSCs in different bacteria evolved from a common ancestor that is not shared by other ion channels. Therefore, the conclusion is that bacterial VGSCs are the most probable channels to be affected by the antagonism. However, there is still no discrete evidence on the best channel among the bacterial VGSC cluster.

Third iteration:
In the next cycle, since we have now clarified only VGSCs are related to each other in phylogenetic relationships, only VGSCs are sampled for phylogeny. 22 proteins including negative control and the new inclusion of VGSCs from archaea are randomly chosen from NCBI database. Archaea VGSCs are chosen because archaea often include proteins that can bridge the relationship between prokaryotes and eukaryotes. In the results in Figure A5, the phylogeny shows that bacillus haloduran VGSCs reside on the eukaryotic cluster.

Figure A5: small_aligned_5_JTT.pdf

We then excluded some distant VGSCs and only isolated those proteins in the same cluster as bacillus halodurans VGSCs as shown in Figure A6.

Figure A6: small_aligned_4_CLUSTALW_JTT.pdf

In the final iteration of the cluster, the cluster showed consistency in the inclusion of Bacillus Halodurans in the eukaryotic VGSCs cluster. So our preliminary choice was Bacillus Halodurans.

However, there are multiple sources pointing out that Bacillus Halodurans is incapable of responding to PSTs, therefore we looked for E coli’s response to PSTs instead, which is demonstrated above.

Below appendix contains equations, table of parameters, wet lab raw data.

References:

1. O. Trott, A. J. Olson, AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading, Journal of Computational Chemistry 31 (2010) 455-461

2. Yang Liu, Maximilian Grimm, et al . CB-Dock: a web server for cavity detection-guided protein–ligand blind docking. Acta Pharmacologica Sinica, 2019.

3. van den Berg, J., Galbiati, H., Rasmussen, A. et al. On the mobility, membrane location and functionality of mechanosensitive channels in Escherichia coli. Sci Rep 6, 32709 (2016).

4. Kenney, Linda J., and Ganesh S. Anand. "EnvZ/OmpR two-component signaling: an archetype system that can function noncanonically." EcoSal Plus 9.1 (2020).

5. Wang, Loo Chien et al. “The inner membrane histidine kinase EnvZ senses osmolality via helix-coil transitions in the cytoplasm.” The EMBO journal vol. 31,11 (2012): 2648-59. doi:10.1038/emboj.2012.99

6. Srividhya, K., Krishnaswamy, S. A simulation model of Escherichia coli osmoregulatory switch using E-CELL system. BMC Microbiol 4, 44 (2004).

7. Pomati, Francesco et al. “Effects of saxitoxin (STX) and veratridine on bacterial Na+ -K+ fluxes: a prokaryote-based STX bioassay.” Applied and environmental microbiology vol. 69,12 (2003): 7371-6.

8. Igoshin, Oleg A et al. “Hysteretic and graded responses in bacterial two-component signal transduction.” Molecular microbiology vol. 68,5 (2008): 1196-215.

9. Barbieri, Christopher M et al. “Comprehensive analysis of OmpR phosphorylation, dimerization, and DNA binding supports a canonical model for activation.” Journal of molecular biology vol. 425,10 (2013): 1612-26.

10. Yoshida, Takeshi, et al. "Transcription regulation of ompF and ompC by a single transcription factor, OmpR." Journal of Biological Chemistry 281.25 (2006): 17114-17123.

11. Yusuf, Rahmi. Synecdoche, Two-Component Systems: E. coli EnvZ/OmpR as. Diss. University of Portsmouth, 2019.

12. Russo, Frank D., and Thomas J. Silhavy. "EnvZ controls the concentration of phosphorylated OmpR to mediate osmoregulation of the porin genes." Journal of molecular biology 222.3 (1991): 567-580.

13. Tigges, M., Marquez-Lago, T., Stelling, J. et al. A tunable synthetic mammalian oscillator. Nature 457, 309–312 (2009).

14. ​​Whitaker, Weston R., et al. "Engineering robust control of two-component system phosphotransfer using modular scaffolds." Proceedings of the National Academy of Sciences 109.44 (2012): 18090-18095.

15. Hsiao, Victoria, et al. "Design and implementation of a biomolecular concentration tracker." ACS synthetic biology 4.2 (2015): 150-161.

16. Saltelli, A. (2002), Sensitivity Analysis for Importance Assessment. Risk Analysis, 22: 579-590.

17. Saltelli, Andrea, et al. Global sensitivity analysis: the primer. John Wiley & Sons, 2008.

18. Zhou X., Lin H., Lin H. (2008) Global Sensitivity Analysis. In: Shekhar S., Xiong H. (eds) Encyclopedia of GIS. Springer, Boston, MA.

19. Zhou X., Lin H. (2008) Local Sensitivity Analysis. In: Shekhar S., Xiong H. (eds) Encyclopedia of GIS. Springer, Boston, MA.

20. Saltelli, Andrea, et al. Sensitivity analysis in practice: a guide to assessing scientific models. Vol. 1. New York: Wiley, 2004.

21. Zhang, X.-Y., Trame, M., Lesko, L. and Schmidt, S. (2015), Sobol Sensitivity Analysis: A Tool to Guide the Development and Evaluation of Systems Pharmacology Models. CPT Pharmacometrics Syst. Pharmacol., 4: 69-79.

22. Sudhir Kumar, Glen Stecher, Michael Li, Christina Knyaz, and Koichiro Tamura (2018) MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms. Molecular Biology and Evolution 35:1547-1549

23. Zakon, Harold H. “Adaptive evolution of voltage-gated sodium channels: the first 800 million years.” Proceedings of the National Academy of Sciences of the United States of America vol. 109 Suppl 1,Suppl 1 (2012): 10619-25.

PROJECT

• Description
• Design
• Engineering
• Proof of Concept
• Implementation
• Safety
• Contribution
• Notebook

WETLAB

• Experiments
• Results

DRYLAB

• Model

REACH

• Human Practice
• Collaboration
• Communication
• Education
• Sustainable

PARTS

• Basic Parts
• Composite Parts

TEAM

• Team
• Attributions