Protein Modelling

Background

Spiders have developed large complex repertoires of peptides in their venoms with various pharmacological and agricultural applications for more than 300 million years (Saez & Herzig, 2019). These peptides are made up of characteristic three-dimensional folds that give them stability and resistance (King, 2019). In this way, spider venom peptides potently and selectively modulate a wide variety of targets, such as ion channels and signaling pathways that are involved in important physiological processes. In addition to the above, due to the multiple disulfide bonds they present, these peptides have remarkable stability within the intestine of insects (King & Hardy, 2013).

Spiders use these peptides to incapacitate prey by modulating the activity of their voltage-gated sodium and calcium channels, which can be classified as pore blockers, which physically obstruct the channel pore, or as modifiers of the channel gate, since they alter the entrance of the channel through interactions with one or more specific domains of the same (Shen et al., 2018).

Selected Peptides

For spider venom peptides selection, an extensive bibliographic research of these was carried out and a database was established with many peptides from different species of spiders. In the first place, the peptides were selected according to their selectivity in invertebrates, and those that showed affinity in vertebrate ion channels were discarded, likewise, an important parameter within the selection of peptides was the insecticidal effect that they exerted on the insect order Coleoptera.

Finally, we reviewed their mechanism of action and the number of disulfide bonds they contained, since these have a great influence on conferring biological stability to the peptide. Based on the specifications mentioned above, the selected peptides were Omega-Hexatoxin-Hv2a (AcTx-Hv2a) and U1-Theraphotoxin-Sp1a (OAIP-1).

AcTx-Hv2a from the Australian spider Hadronyche versuta and OAIP-1 from the Australian tarantula Selenotypus plumipes are two insecticidal peptides that target voltage-gated sodium and calcium channels of insects belonging to the orders Coleoptera, Diptera, Dictyoptera, Orthoptera, and Lepidoptera. Both peptides form a cysteine knot, resulting in three disulfide bonds, providing the peptides with biological stability that contributes to their oral activity and resistance to proteases (Hardy et al., 2013).

OMEGA-HEXATOXIN-HV2A

ω-Hexatoxin-Hv2a (ω-HXTX-Hv2a or AcTx-Hv2) is a specific and potent blocker of voltage-gated calcium channels. It is a 45 amino acid peptide containing a compact disulfide-rich core, including a highly hydrophobic C-terminal tail that is essential for interaction with insect calcium channels, the three disulfide bridges in the hydrophobic core form a binding motif inhibitory cysteine (Wang et al., 2001).

Figure 1. Structure of ω-Hexatoxin-Hv2a. The ammino acids marked in yellow are the ones corresponding to the disulfide bonds.

U1-THERAPHOTOXIN-SP1A

U1-TRTX-Sp1a or OAIP-1 is an insecticidal peptide of 34 amino acids, its three-dimensional structure, like the previous peptide, has six cysteines that form an inhibitory knot formed by three disulfide bonds, providing the peptide an exceptional biological stability (Hardy et al., 2013).

Figure 2. Structure of U1-TRTX-Sp1a. The ammino acids marked in yellow are the ones corresponding to the disulfide bonds.

FUSION PROTEIN

One of the main challenges we encountered using spider venom peptides was their low oral efficiency. These proteins evolved to be injected directly into their prey by the spider. Therefore, they do not possess the mechanisms for oral absorption, being less effective when administered orally. We overcame this challenge by designing a fusion protein, merging each insecticidal peptide with a carrier. Fusion proteins built with a lectin from the Galanthus nivalis agglutinin (GNA) plant significantly improve oral efficiency thanks to its ability to cross through the intestinal epithelium, delivering the peptides to the hemolymph of the insect. Additionally, with the GNA lectin, we improved protein stability and resistance against proteolysis (Powell et al., 2020).

Figure 3. Structure of the lectin from the plant Galanthus nivalis agglutinin.

Figure 4. Mechanism of action of GNA to transport peptides to the insect hemolymph. In the digestive tract, GNA binds to glycoproteins of the intestinal epithelium, crosses the epithelial cell by endocytosis and exocytosis and is transported to the hemolymph.

The fusion proteins assembled were the AcTx-Hv2a/GNA and OAIP-1/GNA. Each peptide was attached to the N-terminal end of GNA with a three-alanine linker, to allow maximum flexibility (Crasto & Feng, 2000). The alanine linker was constructed as reported by Pyati and collaborators, 2014.

Figure 5. Fusion proteins of GNA (dark blue) with a) AcTx-Hv2a and b) OAIP-1 (light blue) via a three-alanine linker (red)

To fuse the AcTx-Hv2a and OAIP-1 peptides with a GNA via the three-alanine linker we used the methodology that follows:

We used the program MOE (Molecular Operating Environment to visualize the sequence of the proteins and create the sequence of the three-alanine linker. We neared the C-terminal ends of the peptides and the N-terminal end of the triple alanine until the corresponding ends were close enough to form a peptide bond. To conduct the bond, we deleted an oxygen atom from the C-terminal of the peptides and a hydrogen atom from the N-terminal of the linker. Subsequently, we were required to change the N-terminal of the linker to a sp2 hybridization and a neutral charge. Once the conditions of the corresponding terminal ends were adequate, we used the molecular editing tool of the software to create the peptide bond between the proteins and the linker. Then, we optimized the bond angles and distances of the binding region with a molecular minimization process.

Finally, to complete the fusion protein, we repeated the methodology to merge the C-terminal of the linker with the N-terminal of the lectin (iGEM UAM, 2021)

Throughout the procedure, we were aided by Dr. Jaqueline Padilla Zúñiga from the Universidad Autónoma Metropolitana.

MOLECULAR DYNAMICS

Molecular dynamics simulations allow modeling the movement of atoms based on Newton's laws of motion, making it an essential tool for drug discovery. One of its many applications is to predict conformation and protein folding (De Vivo et al., 2016).

Previously, we changed the third structure of the proteins when we fused each of them with the GNA carrier. Therefore, we were required to model their new tertiary structure. We used the computational software NAMD and the clusters from the institutional supercomputer YOLTLA to run a two hundred nanoseconds simulation of the solvated fusion proteins at twenty degrees Celcius and a neutral pH.

Figure 6. A) Fusion protein AcTx-Hv2a/GNA pre molecular dynamics analysis B) Fusion protein AcTx-Hv2a /GNA post molecular dynamics analysis C) Fusion protein OAIP-1/GNA pre molecular dynamics analysis D) Fusion protein OAIP-1/GNA post molecular dynamics analysis.

As shown in the pictures above, we can observe that both fusion proteins tend to fold into a globular configuration. Said globularity is more prominent in the OAIP-1/GNA protein than the AcTx-Hv2a/GNA one. We can also notice that the configuration change comes from the folding of the alanine linker and the conjunct amino acid residues from the lectin and the corresponding insecticidal peptide.

MOLECULAR DOCKING

Once we had the model of each protein, we carried out molecular docking analysis between them and their respective targets, the sodium and calcium ion channels. The objective was to validate their action mechanism and hence corroborate their insecticidal effect on the targeted pests. Additionally, it permitted us to verify that the new configurations don't interfere with their insecticidal effect.

We performed the molecular dockings via the server HDock and analyzed each fusion protein with its corresponding target. The protein carrying the AcTx-Hv2a peptide was studied with a voltage-gated calcium channel, and the one with the OAIP-1 peptide with a voltage-gated sodium channel.

Due to the lack of properly characterized coleopteran ion channels, we used the sequence of the ones from the insect Drosophila melanogaster (fruit fly). Said channels have a homology of 89% (Smith et al., 2013) and 74% (King et al., 2008) with the sodium and calcium channels of the coleopteran insect Tribolium castaneum, respectively, providing suitable substitutes for the study.

We evaluated the molecular docking results with the specialized software Discovery Studio. This program allowed us to visualize and analyze the interactions between the fusion proteins and the ion channels and, subsequently, to compare the results with the ones reported in the literature.

RESULTS & ANALYSIS

Ω-HEXATOXIN-HV2A/GNA—CAV CHANNEL

As previously said, the ω-Hexatoxin-Hv2a is a potent pore blocker of the voltage-sensitive calcium channels. The action mechanism of this highly specific peptide is not yet fully comprehended. Nonetheless, there are putative mechanisms reported in the literature. The essence of the action mechanism is for the disulfide-bond-rich zone of the peptide to interact with the extracellular residues of the ion channel, interfering with the flow of ions (Wang et al., 2001).

Figure 7 .A) Side view and B)Top view of the interaction of the AcTx-Hv2a fusion protein (dark and light blue) and the calcium channel (violet/purple).

As seen in the previous images, the AcTx-Hv2a sites at the entrance of the calcium channel, obstructing the passage of ions. Distinct from the literature, the blocking of the pore is conducted by the insecticidal peptide and the carrier. Additionally, we can appreciate that the fusion protein does not entirely block the flow of calcium. However, it reduces enough the passage of ions to cause an inhibitory effect.

OAIP-1/GNA—NAV CHANNEL

The OAIP-1 is an inhibitory peptide that targets the sodium ion channels. This action method is known to cause a shift in the polarization of the ion channel, modifying the passage of ions (Smith et al., 2013). More specifically, the insecticidal peptide targets the site 4 of the voltage-sensing domain II (King, 2019). By engaging with experts, we discovered that the site 4 of the VSD has a series of positively charged amino acids (Arg) which are the essence of the sensor domain. Furthermore, any change in at least one of these residues causes the depolarization of the channel resulting in the death of the insect.

The molecular dockings simulation with the OAIP-1/GNA protein and the NaV channel demonstrates that the fusion protein conserves its insecticidal effect as the lectin does not interfere with the action method of the OAIP-1.

Figure 8.Interaction between the fusion protein (dark and light blue) and the VSD II (pink) of the sodium channel

Figure 9. Zoom into the attractive charged interaction of the ARG 909 of the NaV channel with the ASP 12 of the OAIP/GNA (yellow).

As shown, the OAIP-1/GNA interacts with the VSD II of the sodium channel, as reported in the literature. More specifically, the asparagine 12 of the fusion protein engages with the arginine 909 of the sodium channel via an attractive charged interaction.

ENGAGING WITH STAKEHOLDERS

PROTEIN MODELING EXPERT

All this work would not have been possible without the help and advice of experts in the area of protein modeling.

One of the most important elements in our project was the design of our fusion protein, to provide oral efficacy to spider venom peptides. We owe that construction to Dr. Jaqueline Padilla Zúñiga, a researcher at the Department of Chemistry of the UAM Iztapalapa Unit, with her help we were able to design and assemble a fusion protein for each of the insecticidal peptides. This construction was carried out in the specialized MOE software, whose procedure is described above.

Figure 10. Dr. Jaqueline Padilla Zoom meeting

ION CHANNEL EXPERT

Just as it is important to know in depth the peptides with which we are working, it is essential to know the mechanism of action of these peptides with their receptor, in this case, the ion channels of insects. Each peptide has a different action mechanism, which is why we engaged in talks with Dr. Erika Puente Guzmán, an expert in the functioning of ion channels at the UAM Cuajimalpa Unit. With her mentoring, we were able to understand the action mechanism of our peptides and how they exert their insecticidal function on the ion channels of insects.

Figure 11. Dr. Erika Puente Zoom meeting

References

  • Crasto, C. J., & Feng, J. A. (2000, May). LINKER: a program to generate linker sequences for fusion proteins. Protein Engineering, Design and Selection, 13(5), 309–312. https://doi.org/10.1093/protein/13.5.309

  • De Vivo, M., Masetti, M., Bottegoni, G., & Cavalli, A. (2016, February 8). Role of Molecular Dynamics and Related Methods in Drug Discovery. Journal of Medicinal Chemistry, 59(9), 4035–4061. https://doi.org/10.1021/acs.jmedchem.5b01684

  • King, G. F., & Hardy, M. C. (2013). Spider-venom peptides: structure, pharmacology, and potential for control of insect pests. Annual review of entomology, 58, 475-496

  • King, G. F., Escoubas, P., & Nicholson, G. M. (2008, March 5). Peptide toxins that selectively target insect NaVand CaVchannels. Channels, 2(2), 100–116. https://doi.org/10.4161/chan.2.2.6022

  • King, G. F. (2019). Tying pest insects in knots: the deployment of spider-venom-derived knottins as bioinsecticides. Pest Management Science, 75(9), 2437–2445. https://doi.org/10.1002/ps.5452

  • Pyati, P., Fitches, E., & Gatehouse, J. A. (2014, August 1). Optimising expression of the recombinant fusion protein biopesticide ω-hexatoxin-Hv1a/GNA in Pichia pastoris: sequence modifications and a simple method for the generation of multi-copy strains. Journal of Industrial Microbiology and Biotechnology, 41(8), 1237–1247. https://doi.org/10.1007/s10295-014-1466-8

  • Saez, N. J., & Herzig, V. (2019). Versatile spider venom peptides and their medical and agricultural applications. Toxicon, 158, 109–126. https://doi.org/10.1016/j.toxicon.2018.11.298

  • Shen, H., Li, Z., Jiang, Y., Pan, X., Wu, J., Cristofori-Armstrong, B., ... & Yan, N. (2018). Structural basis for the modulation of voltage-gated sodium channels by animal toxins. Science, 362(6412)

  • Smith, J. J., Herzig, V., King, G. F., & Alewood, P. F. (2013, March 23). The insecticidal potential of venom peptides. Cellular and Molecular Life Sciences, 70(19), 3665–3693. https://doi.org/10.1007/s00018-013-1315-3

  • Wang, X. H., Connor, M., Wilson, D., Wilson, H. I., Nicholson, G. M., Smith, R., Shaw, D., Mackay, J. P., Alewood, P. F., Christie, M. J., & King, G. F. (2001, October). Discovery and Structure of a Potent and Highly Specific Blocker of Insect Calcium Channels. Journal of Biological Chemistry, 276(43), 40306–40312. https://doi.org/10.1074/jbc.m105206200

Mathematical Modelling

Introduction

Synthetic biology is an engineering discipline that builds on our mechanistic understanding of molecular biology to programming microbes to carry out new functions. Such predictable modification of a cell requires modelling and experimental techniques to work together. The modelling component of synthetic biology allows one to design biological circuits and analyze their expected behavior. The synthetic merges models with evident systems by providing quantitative data and sets of available biological parts that can be employed to construct circuits. Sufficient progress has been made in the combined use of modelling and experimental methods, which reinforces the idea of being able to use engineered microbes as a technological platform (Chandran et. al., 2008).

In other words, for iGEM teams, the mathematical model is a useful tool that can make the wet lab easier to improve the planned experiment. It consists of trying to synthesize the nature of a phenomenon from computational tools and does not try to demonstrate the phenomenon (iGEM Foundation, 2022). In this way, we can do a simple or idealized description of our system, situation, or process, often in mathematical terms, that is put forward as a basis for theoretical or empirical understanding or for calculations, predictions, etc. (OED Online, 2022).

From mathematical equations we can build a mathematical model that is a set of differential equations. Then one can go from biochemical reactions to equations that indicate how the different concentrations of the species in the system change as a function of time. There are even more reasons why the mathematical model is relevant for synthetic biology (iGEM Foundation, 2022):

  • Prove that an experimental design is realizable

  • It guides us in the selection of parts and assembly

  • Predict expected results before experiments in the lab

  • Give us feedback so we can improve the design, redesign if is the case

  • Run experiments under conditions that may be dangerous or costly in the real system

  • It shows if it is necessary to reduce the risk inherent in decision making

Background

Our model represents the interaction mechanism between the genetic circuit elements and all the proteins that interact with the latter over time. Since we could not transform the Pichia correctly with our fusion protein to do experiments, we tried to make useful predictions from a model, using the parameters accurately estimated. Mechanistic models were established on the fundamental physical and chemical laws. Including parameters that carry physical, chemical, or biological principles. For this, less experimental data are necessary to determine parameters and calibrate the model (Zheng & Ganesh, 2010)

The parameters that define our system can be found in literature, experimental, or screening. For literature, iGEM teams have previously registered constants that can help us, so we must find at least a range of values; experiments refer to when parameters are relevant to measure in the laboratory; screening is required when we do not found any data in the literature. Then we proceed to perform a screening of parameters (iGEM Foundation, 2022). In this case, we implemented the screening parameters.

Last year we performed a sensitivity analysis and implemented a genomic scale metabolic model; the mathematical model was then adjusted and helped us simulate the behavior of systems to optimize the production of Spidicide-CX. However, the results were not specific to our system. Then, we took a rank to predict relevant parameters for our system. For this, we based on previous iGEM teams with similar parameters we obtained. Those are iGEM Colombia’s Team (2012, 2013); iGEM iBowu-China's Team (2019) .

After speaking with our advisers, we knew that was necessary to find a specific value of just one or two parameters to make it simple. To create a prediction of the model closer to what we could obtain by having our transformed cell ready and performing the respective experiments. This value would be looked at in papers where the aim would be like what we were looking. We found only one parameter in the literature that, after performing the screening, we knew that it was one of the parameters sensitive to the model.; we could not have the values we wanted experimentally, but we had another solution to solve this situation (hardware and software).

Mechanistic model

We can define a mechanistic mathematical model as the mathematical description of the elements forming a system, their mutual interactions, and their interaction with the environment. Such models are used in technical systems to enable the extrapolation of systems' behavior relying on the mathematically described features of elements and mechanisms of their interaction (OED Online, 2022).

Modelling Description

For our model we employ differential equations to express the mass balances of the system, the representation of the molecular mechanisms uses what is known as the law of mass action and Hill's equation

Where:

  • 𝛽0 is the basal expression (if it is 0 there are no activators)

  • A is the number of activators

  • n is Hill's coefficient

  • d is the association constant

These equations complement each other and use the experimental parameters obtained in literature to represent the relationships of the compounds involved.

Genetic circuit to produce OAIP-1 and AcTx-Hv2 in Pichia pastoris

In this case, we focus on a scenario where the main objective is the production of the fused protein (FP) OAIP-1 and AcTx-Hv2. As we know, they are produced separately but have the same genetic circuit. The two proteins are presented as PF but remember the last mentioned.

Figure 1: Genetic circuit illustrated

Enzymatic reactions

First, of all we would represent the enzymatic reactions that are happening inside the cell:

Figure 2: Enzymatic reactions

Differential equations

From those enzymatic reactions, differential equations were proposed. Differential equations express the mass balances of the system, and the representation of the molecular mechanisms uses what is known as the law of mass action and Hill's equation. These equations complement each other and use the experimental parameters obtained in the literature to represent the relationships of the compounds involved. With the help of this equation, it is possible to describe the behavior of each element of the events that are happening inside the cell.

All the equations represented below are based on the law of mass conservation:

Accumulation = Input - Output + Generation - Consumption

Below are the equations and their respective description of what happens in our model for the fused protein production.

Table 1. Differential equations

Symbol

Description

Equation meaning

Equation

MetHOout Methanol as an inducer of the circuit. Its concentration only depends on the diffusion rate. Out-of-cell concentration depending on the diffusion rate to enter or leave the cell.
MetHOin The amount of methanol inside the cell depends on its degradation rate, diffusion rate, and binding to transcription factors. Circuit and amount of methanol.
mRNA The concentration of mRNA depends on the basal production of this molecule. Degradation rate of mRNA
PFGFP Production of the fusion protein and GFP mediated by translation and diffusion out-of-cell rate. Rate of degradation within the cell
PFout Concentration of fusion protein outside the cell depending on the diffusion rate of the cell. Rate of degradation outside.

Parameters

Parameters are the key and limiting step of model construction; this quantity is fixed in a particular case and can vary in different situations. Without them, we cannot proceed to demonstrate something we want to prove. If an initial model exhibits significant departures from experimental data, further experiments may need to be accomplished to refine the parameter values. The process is iteratively repeated until the construction of a satisfactory model.

There are many ways we can obtain the parameters; it could be from literature or even estimated manually. Frequently, the estimation of the parameters is obtained from experimental observations. In our case, that was not possible. We are still experimenting with our chassis. Consequently, we have not yet acquired a product of interest to carry out such tests. Amongst this parameter estimation, a computationally intensive but powerful and robust method is used.

Given the lack of data, we take a range of values to test which parameters are relevant to our project. We get these parameters from other iGEM teams since these parameters have not been found in the literature. Below are the parameters that we believe are the most relevant. Those that do not show significant sensitivity are neglected.

Table 2. Parameters

Estimated range of values of the parameters involved in the mathematical model (iGEM Colombia’s Team, 2012, 2013; iGEM iBowu-China's Team, 2019; Valencia UPV iGEM, 2018)

Description of sensitivity analysis

In a numerical model, the Sensitivity Analysis (SA) is a method that measures how the impact of uncertainties of one or more input variables can lead to ambiguity on the output variables. This analysis is appropriate because it improves the prediction of the model. Reduces it by studying qualitatively and quantitatively the model response to a change in input variables or by understanding the phenomenon studied by the analysis of interactions between variables (Wexler P., 2014).

In the development of our project, the production of a fusion protein implies more characteristics and behaviors than previously known. In addition, the lack of information about it in the literature. Even though determining specific parameters by experimenting was hard to execute. That is why we change the approach of the modelling part.

At the first step of design, a sensitivity analysis was performed using the specialized MATLAB software to determine the importance of each parameter in our fusion protein production. The results allow us to focus on the parameters that will be needed to be experimentally measured in the system to enhance in a second iteration the mathematical model. Those will help us with simulations and predictions, previously made in silico, geared toward biopesticide production.

The range for each parameter was taken from the literature and set to a specific value. Hence, each parameter established in the model has performed the sensitivity analysis. Afterward, the value of the parameters analyzed changed according to the various signal. The concentration of fusion protein inside the cell is the response signal, and the one that reaches the steady state is what the analysis determined.

Obtaining experimental data

Our models need an adjustment to improve them. In this way, we need to obtain the parameters which describe our system. But the specialized literature offers few reported parameters that are useful to us because our chassis is uncommon.

We can get the parameters which we need by experimentation. The particularities of our strain of P. pastoris would require experimentation, in which the kinetic parameters are studied to adequately fit our mathematical model. However, the equipment necessary for such a task was not easily accessible to us and to solve this problem we decided to create hardware capable of performing some of these measurements.

We didn't create a revolutionary device; we made a homemade spectrometer. The Aachen team from Germany made a very similar device in 2014. We proposed to develop a visual interface to be more user-friendly.

Our device seeks to measure optical density and fluorescence using cheap and easily accessible materials. We rely on the Arduino platform, in conjunction with the python programming language to build the necessary hardware and software. The use and design of hardware and software will be covered in more depth in this section.

Our objectives for this device are the followings:

  • Make a homemade spectrophotometer that can monitor the fluorescence and optical density of a bioreactor in real-time.
  • Make a script in Python, which allows serial communication with Arduino. Transmission of signals to LEDs and reading of sensors.
  • Make a circuit with an Arduino, which allows serial communication with Python, which contains 2 LED outputs and the reading of two light sensors.
  • In the Python code, we must include a graphical interface which has the following options: Connect to the Arduino circuit, create an optical density measurement mode, create a fluorescence measurement mode, save the measurements to a CVS file, and display the data in a graph in real-time.
  • Simulate the connection in real-time with a replica of the Arduino circuit in Proteus.

You can obtain the protocol of our parameter clicking the following

Parameter protocol

Note: we decided to calibrate our device using E. coli instead of P. pastoris because E. coli, which also produces GFP, was already transformed when the device was developed, thus, the protocol was elaborated in a hypothetical aspect.

Metabolic model

Genome-scale metabolic model (GEM) of Pichia pastoris

We implemented a genomic scale metabolic model that incorporates the production of our biopesticide to design growth strategies for our microorganism at the bioreactor level. Based on this, we adjust the mathematical models that help us simulate the behavior of systems to optimize the production of the Spidicide-CX.

GEM reconstruction has been established as one of the main modelling approaches for systems-level metabolic studies. GEM allows metabolic forecasting flux values for an entire set of metabolic reactions using optimization techniques, such as flux balance analysis (FBA), which uses linear programming (Gu C. et al., 2019). GEMs are a network-based tool that collects all known metabolic information of a biological system, including the genes, enzymes, reactions, associated gene-protein-reaction (GPR) rules, and metabolites. These metabolic networks provide quantitative predictions related to growth or cellular fitness based on GPR relationships.

Consequently, from the mechanistic model, it is possible to study a potential behavior of the design according to changes in the protein due to the concentration of certain factors such as transcribers, and inducers, among others. Information about the Pichia that performs the product transformation, can be added to optimize the production conditions of the peptide of interest. In this way, the metabolic model can be switched to the genomic scale and finally predict the rate of protein synthesis, before experimenting with the wet lab.

To simulate the production of OAIP-1 and AcTx-Hv2a, we followed the methodology of Sohn et al., and the Genome-Scale Metabolic Model of Pichia pastoris, taken from the BioModels repository (Tomàs-Gamisans, 2016) (Sohn S. et al., 2010). Three products were taken as a reference to represent the production mechanism of each protein: their respective DNA, mRNA, and amino acid sequence.

Table 3. Metabolic reactions representing the production of fusion proteins OAIP-1 and AcTx-Hv2a

With the help of the MATLAB Cobra Toolbox tool and based on these equations, simulations were developed by selecting the maximization of the biomass reaction to obtain an optimization.

Through this analysis, the optimal rate of carbon consumption for each type of protein produced was resolved. Next, the impact of methanol consumption was analyzed when the System works correctly to produce the fusion protein. The results could be represented with a heat map where the horizontal axis indicates the carbon consumption, the vertical axis the protein production rate and the specific growth rate through the contour lines.

Although the analysis of the rate of carbon consumed was not enough, it was intended to analyze the oxygen consumption from the optimum values of the inducer. This to demonstrate what the literature indicates, P. pastoris in a hypoproteins. Therefore, the production of fusion proteins was determined, as well as the growth rate at different oxygen concentrations.

Results

Mechanistic model

Now, the results obtained from the sensitivity of the mechanistic model are presented below. It should be noted that the results were divided into three levels of sensitivity: none, low and high. None refers to the fact that there is no slope on the graph, it is a constant value. Low is considered when the slope of the graph is minimal and high when the slope is very noticeable and there is a drastic change in the concentration of the fusion protein.

Table 4. Results of the mechanistic model

A high sensitivity is equal to say that a minimal change in this parameter will drastically affect the intracellular concentration of fusion proteins. This study allows us to deduce which parameters are important to measure in the laboratory and consequently obtain better predictions from our model.

Within the parameters analyzed, there is only high sensitivity in mRNA degradation rate, translation rate, FP1/FP2 diffusion rate and FP1/FP2 degradation rate.

According to the parameters that had the greatest effect on the mechanistic model, we performed a simple optimization analysis within the previously established ranges. Due to the lack of experimental information, an estimation approach for the parameters with the greatest effect can be performed to find the values at which we can find a higher concentration of the fusion protein. In the future we will be able to show whether our estimates are correct by fitting the model to experimental data.

The following graph shows the results of the simulation with the sensitive parameters that estimate the highest concentration of fusion protein in the shortest time. This simulation, like the previous ones, was done in Matlab using the ode45 function to solve the differential equations.

Figure 3: Model prediction using the best sensible parameter values

Over time we can observe a decrease in the concentration of methanol outside the cell, due to its consumption as an inducer in the production of the fusion protein. This simulation presents a maximum peak and the highest production rate of the PF in the first 15 minutes of the process, reaching a semi-steady state after 35 minutes. This simulation represents the best capacity of our circuit using the best sensitive parameters found in the sensitivity analysis, whose specific values are: mRNA degradation rate 0.3 min-1, translation rate 0.7 min-1 and PF degradation rate 0.05 min-1.

Metabolic model

The simulations results using the previous described GEM on P. pastoris showed the carbon flux distributions for both fusion proteins.

A heatmap data visualization is a method of graphically representing numerical data where the value of each data point is indicated using colors. For OAIP-1 and AcTx-Hv2a production, the same behavior is observed. The heatmap allowed us to obtain a range of values of methanol and oxygen consumption rates.

For both fusion proteins, the rate of P. pastoris is more important when the consumption of methanol is high and protein production rate is low. This is due because the protein production is mutually exclusive with the specific growth rate. A curious fact is that we could achieve a considerable number of fusion proteins without having biomass. One thing we must consider is that at high concentrations of methanol, it can damage the cell, so we decided to consider a rate of methanol consumption rate value of 7 mmol/gDCW/h.

An optimum performance is present at a low oxygen uptake rate for both cases of fusion protein production. In this situation, the highest values of the protein or biomass production rate can be reached.

We can then conclude that the best conditions in which there is an optimal production of protein and biomass is at a rate of methanol consumption of 7 mmol/gDCW/h and oxygen consumption rate of 2 mmol/gDCW/h.

OAIP-1

AcTx-Hv2

Methanol uptake rate

Oxigen uptake rate

Figure 4. Heat map results of flux variability analysis of OAIP-1 and AcTx-Hv2a

Finally, we can say that this analysis of the metabolic model together with the mechanistic model, achieve a good prediction to reach an optimal value of fusion protein production. Therefore, the values obtained by mathematical modelling will facilitate laboratory experiments.

Future Vision

In its broadest sense, mathematical modelling is the process of a problem solving by the mathematical expression of real-life event or a problem. The study we have done have already done a great prediction. However, laboratory measurements of the parameters that were considered sensitive for protein production have yet to be taken. This is one the priority will be to know the parameters of the model by designing the required experiments for next year.

It is also planned to carry out measurements with the transformed Pichia sensor. Recall that protein production is equivalent to the presence of fluorescence. Therefore, we could verify the results of the metabolic modelling, that with high concentrations of methanol and low concentrations of oxygen there will be greater fluorescence. Likewise, it would be interesting to verify how much protein is produced when a methanol shock is administered to the cell medium.

The protocol that will be necessary to follow to verify these assumptions is presented below.

Figure 5. Protocol for obtaining model parameters

REFERENCES

  • Chandran, D., Copeland, W. B., Sleight, S. C., & Sauro, H. M. (2008). Mathematical modelling and synthetic biology. Drug discovery today. Disease models, 5(4), 299–309. https://doi.org/10.1016/j.ddmod.2009.07.002

  • iGEM Foundation [iGEMFoundation]. (2022, july 14). Modelling for SynBio: from ODEs to Gene expression. Youtube. https://www.youtube.com/watch?v=ph5iYWwXsPw

  • "model, n. and adj." OED Online. Oxford University Press, September 2022. Web. 4 October 2022.

  • Team:Colombia Uniandes/Parameters - 2013.igem.org. (s. f.). Recuperado 4 de octubre de 2022, de https://2013.igem.org/Team:Colombia_Uniandes/Parameters

  • Stalidzans, E., Zanin, M., Tieri, P., Castiglione, F., Polster, A., Scheiner, S., Pahle, J., Stres, B., List, M., Baumbach, J., Lautizi, M., Van Steen, K. & Schmidt, H. H. (2020, 1 mayo). Mechanistic Modelling and Multiscale Applications for Precision Medicine: Theory and Practice. Network and Systems Medicine, 3(1), 36-56. https://doi.org/10.1089/nsm.2020.0002

  • Randers-Eichhorn, L., Albano, C. R., Sipior, J., Bentley, W. E., & Rao, G. (1997). On-Line Green Fluorescent Protein Sensor with LED Excitation. In John Wiley & Sons, Inc. Biotech-nol Bioeng (Vol. 55).

  • Team:Alemania Recuperado 6 de octubre de 2022, de https://2014.igem.org/Team:Aachen/OD/F_device

  • Wexler, P. (2014). Encyclopedia of Toxicology. Elsevier Gezondheidszorg.

  • Gu, C., Kim, G.B., Kim, W.J. et al. Current status and applications of genome-scale metabolic models. Genome Biol 20, 121 (2019). https://doi.org/10.1186/s13059-019-1730-3

  • Passi, A., Tibocha-Bonilla, J. D., Kumar, M., Tec-Campos, D., Zengler, K., & Zuniga, C. (2021). Genome-Scale Metabolic Modelling Enables In-Depth Understanding of Big Data. Metabolites, 12(1), 14. https://doi.org/10.3390/metabo12010014

  • Tomàs-Gamisans2016 - Genome-Scale Metabolic Model of Pichia pastoris (version 2) | BioModels. (s. f.). Recuperado 8 de octubre de 2022, de https://www.ebi.ac.uk/biomodels/MODEL1508040001

  • Sohn, S. B., Graf, A. B., Kim, T. Y., Gasser, B., Maurer, M., Ferrer, P., Mattanovich, D., & Lee, S. Y. (2010). Genome-scale metabolic model of methylotrophic yeast Pichia pastoris and its use for in silico analysis of heterologous protein production. Biotechnology journal, 5(7), 705–715. https://doi.org/10.1002/biot.201000078.

  • Principle and Feature of Various Detection Methods (2) : Hitachi High-Tech GLOBAL. (s. f.). Recuperado 9 de octubre de 2022, de https://www.hitachi-hightech.com/global/products/science/tech/ana/lc/basic/course8.html

  • Myers JA, Curtis BS, Curtis WR, et al. Improving accuracy of cell and chromophore concentration measurements using optical density. BMC Biophys. 2013;6(1):4. doi:10.1186/2046-1682- 6-4.

  • Kiviharju K, Salonen K, Moilanen U, Eerikäinen T. Biomass measurement online: The performance of in situ measurements and software sensors. J Ind Microbiol Biotechnol. 2008;35(7):657-665. doi:10.1007/s10295-008-0346-5.

  • Marinescu, Catalin & Popescu, Roua. (2018). Open-Source bioreactor controller for bacterial protein expression. 10.7287/peerj.preprints.27150v1.