Model
“In order to solve a differential equation you look at it till a solution occurs to you.” - George Pólya
Overview
Methods
Results
Introduction
Conclusion
References
Introduction
Mathematical models in science are often used to support experiments in conjunction with theory. By simulating experiments with models’ new insights can be found, and it can direct new experiments. Consequently, fewer experiments can be done, and it increases the opportunity for useful results. By replacing certain experiments with simulations, advantages are created. Simulations are less dangerous, have a high degree of control (just change one parameter), can test higher extremes than in the laboratory, are cheaper, and observe processes that cannot be easily observed by experimenting, such as the progression of a pathway.
We made a kinetic model about the signaling pathway of our engineered cell to simulate the generalized extracellular molecule sensor (GEMS) receptor and its activation. This model is designed to simulate the behavior of the system in normal and in extreme situations by only changing certain variables. This way, the system’s behavior can be predicted without getting dangerous situations in the laboratory. Furthermore, the model especially supports us in improving and reducing future lab work by predicting optimizing conditions and useful experiments.
The goal of our project, !MPACT, is to improve upon the current unspecific treatments against ANCA-associated vasculitis (AAV) by creating a Modular and Personalized Autoimmune Cell Therapy. First of all, the engineered !MPACT cells should be able to bind antineutrophil cytoplasmic antibodies (ANCAs), disease markers of AAV, and produce a therapeutic protein that can treat AAV. ANCAs bind to the GEMS receptor with as extracellular domain the proteinase 3 (PR3) protein, the target protein of ANCAs. The binding of ANCAs to the receptor activates the JAK/STAT signaling pathway, which in turn leads to the transcription, translation, and secretion of a therapeutic protein. Based on discussions with stakeholders, and literature research, we concluded that interleukin 10 (IL-10) fits the purposes of the therapeutic protein. Therefore, IL-10 became the output of our engineered cells. There were some doubts about using IL-10 as a therapeutic protein by the University Medical Centers of Rotterdam and Utrecht, see Milestone 4. Since the immune system is extremely well-regulated, small changes, such as increasing the IL-10 concentration, can have major consequences on the immune system which might lead to undesired effects. Therefore, it is important that the concentration of IL-10 can be controlled and predicted. Here our model comes into play.
The biochemical reactions in !MPACT cells, from the binding of ANCA to the PR3 protein of the GEMS receptor to IL-10 secretion, are simulated with an ordinary differential equation (ODE) model. The pathway that is simulated can be seen in Figure 1. ODE models consist of a set of differential equations that describe non-linear dynamic systems over time.1 An ODE model can simulate, for example, signaling pathways,2 biochemical reaction networks,3 or gene regulation4 as a function of time. We aim to create a better understanding of the behavior of the engineered GEMS system. With an ODE model, our signaling pathway can be simulated and the behavior of the system, in the form of concentration graphs over time, can be predicted.
The goal of the model is to support future lab work by predicting experimental results, such as the IL-10 concentration, and finding sensitive parameters that can be modified in the laboratory. By modifying the parameters in the model, the system can be optimized to an IL-10 concentration within the therapeutic window. In addition, modifying the sensitive parameter values and initial concentrations can reduce the number of experiments and optimize the conditions of the experiments. This allows the experimental lab work to be used to its full potential. In this way, we want to support our further research and also other researchers using the GEMS system during their lab work. This is achieved by the modularity of our model because a different activation of the pathway and protein secretion can be easily introduced.
With the ODE model of !MPACT, we were able to simulate the engineered GEMS system and after validation of the model with literature we were able to predict the production of IL-10. Also, a multiple parameter sensitivity analysis (MPSA) was performed to explore which parameters are the most sensitive and influence IL-10 production the most. We suspect that changing the downstream parameters will have a greater influence on the IL-10 concentration over time than changing the upstream parameters. Whereas a signal pathway undergoes amplification, a change in the rate of the biochemical reactions at the end of the pathway will have a greater influence on the IL-10 production than a reaction at the beginning. From the results of the MPSA, the optimal values for the parameters are chosen.
The ODE model, with the parameter values of the MPSA, is able to predict dose-response behavior, the maximum limit of response, and the time between binding of ANCAs to our GEMS receptor and the maximum production of IL-10.
Some of the simulations are also performed with the ANCA concentration present in AAV patients. This ANCA concentration is used as ligand concentration when the optimal ratio between ligand and receptor was determined. We suspect that the optimal ratio between ligand and receptor will be approximately 1:100 (ligand: receptor) based on previous results. This ratio can be used in the other simulations. We aimed to use this as a suggestion for the next step in the lab to make our proof of concept easily applicable to the human body.
We also look at how antibodies with different dissociation constants (KD), between antibody and receptor, will influence IL-10 production over time. We expect the concentration of IL-10 to increase as the KD decreases. With these results, the receptor of the cell can be designed, in the future, in such a way that the ANCAs bind with a KD that provides IL-10 secretion within the therapeutic range.
The results predicted with the model can define future experiments with the GEMS platform coupled to the JAK/STAT signaling pathway. The model is easily adjustable, so, when a different receptor or therapeutic protein is used, several parameters and initial concentrations need to be changed while the other parameters are constant.
Conclusion
All the tests to determine whether the model satisfied biologically feasible behavior were successfully accomplished. This demonstrates that the model exhibits biologically feasible behavior. The model achieves equilibrium, which is desirable because almost all biological processes achieve a steady state for their concentrations.
For the MPSA, we expected that KD1 or KBu11 would have the biggest influence on IL-10 secretion because they influence the binding time. KBu11 was the most sensitive, by slowing down this rate the IL-10 production can be increased and the other way around. The laboratory should focus on the duration of the binding between ANCAs and receptors if the secretion of IL-10 needs to be influenced or optimized.
A maximum receptor-ligand ratio of 500 was deduced. Increasing the receptor concentration to an even higher receptor-ligand ratio does not affect the IL-10 concentration and is therefore ineffective. To get a dimerized receptor and pathway activation, the effective receptor-ligand ratio had to be at least 5.
The MPSA results showed that varying the KD1 value at 36 hours had no effect on IL-10 concentration. When the IL-10 concentration was examined at different time points for the range of KD1 values, it was discovered that a picomolar range causes a faster increase in IL-10 concentration and a micromolar range causes a slower increase in IL-10. To summarize, KD1 affects the IL-10 concentration, but the effect is time-dependent.
- Alexandre Megretski. Lecture 2 6.243J Dynamics of Nonlinear Systems, Fall 2003. Dynamics of nonlinear systems (2003).
- Bachmann J, Raue A, Schilling M, Becker V, Timmer J, Klingmüller U. Predictive mathematical models of cancer signalling pathways. In: Journal of Internal Medicine. Vol 271. ; 2012. doi:10.1111/j.1365-2796.2011.02492.x
- Neelamegham S, Liu G. Systems glycobiology: Biochemical reaction networks regulating glycan structure and function. Glycobiology. 2011;21(12). doi:10.1093/glycob/cwr036
- Chai LE, Loh SK, Low ST, Mohamad MS, Deris S, Zakaria Z. A review on the computational approaches for gene regulatory network construction. Comput Biol Med. 2014;48(1). doi:10.1016/j.compbiomed.2014.02.011
Model structure
MPSA
References
Model structure
The Ordinary Differential Equation (ODE) model was created in MATLAB and solved with the ODE15s solver. Our ODE model is based on an existing model of the Generalized Extracellular Molecule Sensor (GEMS) platform with the JAK/STAT signaling pathway made by Jolien Marcelis (our team captain), Stijn van Wijngaarden, Tibo Verburg, and Anton van Galen during their Bachelor End Project. This general model was based on coiled coils as input, the JAK/STAT signaling pathway, and SEAP over time as output. The activation of this GEMS platform is adjusted to the binding of antibodies to proteinase 3 (PR3) and the production and secretion of SEAP is interchanged to interleukin 10 (IL-10).
Description of the simulated system
An overview of the model is presented in Figure 1. Our JAK/STAT3 signaling model is based on the JAK/STAT1 signaling model by Yamada et al. (2003) as the signaling pathways are conserved in cells and STAT1 and STAT3 have the same cellular processes.1-3 In our system, the pathway is activated by anti-neutrophil cytoplasmic antibodies (ANCAs) binding to the GEMS receptor. STAT3 is phosphorylated when bound to the receptor and gets dimerized. Dimerized STAT3 can transfer to the nucleus, where it functions as a transcription factor for IL-10 expression. The IL-10 mRNA is produced and, eventually, the IL-10 protein is formed. When the STAT3 gets dephosphorylated it moves from the nucleus to the cytoplasm where it can bind again to the GEMS receptor.
Receptor activation
The ligand (L) is an anti-neutrophil cytoplasmic antibody (ANCA) and binds to the extracellular domain of the GEMS receptor, which consists of proteinase 3 (PR3) protein bound to an EpoR subunit. JAK is already intracellularly bound to the receptor (RJ).4,5 Binding of ANCA (L) to the receptor subunit (RJ) forms a ligand-receptor complex (RJL), step 1.4,5 After binding the ligand to one receptor subunit, this receptor is able to dimerize with another receptor subunit (forming RJ2L). Rotation of the receptor subunits, caused by dimerization, alters the intracellular architecture.6 Due to this alteration, JAKs are activated by phosphorylation, step 2. This receptor complex is able to phosphorylate several tyrosine residues at the docking sites of the intracellular part of the RJ2L complex, resulting in a phosphorylated dimerized receptor complex called RJ2Lp.4
STAT3 in the cytoplasm (STAT3c) can bind to the docking sites of RJ2Lp that arise from the phosphorylated tyrosines, and form a complex called RJ2Lp-STAT3c, step 4. STAT3c gets phosphorylated by JAK, which causes the dissociation of STAT3cp (phosphorylated STAT3c) from the RJ2Lp, step 5.4 This reaction is reversible, so STAT3cp can again bind to RJ2Lp to form the RJ2Lp_STAT3cp and inhibit the binding of STAT3c to RJ2Lp (Figure 2).
STAT3 cycle
STAT3cp can form a homodimer (STAT3cpd), step 6 of Figure 1, that actively moves to the nucleus. In the nucleus, it is called STAT3npd. Here it can function as a transcription factor, or it can break down into two STAT3np monomers. When STAT3npd is broken down into two monomers. These monomers get dephosphorylated to STAT3n by phosphatases in the nucleus (PPN). After which, the STAT3n monomers return to the cytoplasm to become STAT3c again and restart their cycle (Figure 3). Phosphatases are not only present in the nucleus, but also in the cytoplasm. The phosphatases in the cytoplasm are called PPX, step 12 of Figure 1.1,7
IL-10 production
When the dimer STAT3npd is not dissociated into monomers, it functions as a transcription factor of the pLS13-IL-10 plasmid, step 7 of Figure 1. Hence, the transcription of this plasmid gives mRNA of IL-10 (mRNA_IL10). This mRNA gets translated into an IL-10 protein (IL10er) in the Endoplasmic Reticulum (ER) by ribosomes, in step 9, or the mRNA_IL10 gets degraded, in step 8. The IL-10 protein from the ER moves to the Golgi Complex (IL10g), step 10, before it is secreted in the extracellular environment (IL10ex), step 11.
Model construction
The model consists of twenty-two ordinary differential equations that describe the concentration balance of individual proteins or protein complexes. These ODEs are derived from the chemical reaction equations which describe the pathway. The ODEs are implemented in MATLAB and solved with the ode15s function. The MATLAB scripts can be found on GitHub.
All reactions are represented by mass-action kinetics, except for the transcription and translation reactions, which are based on Hill kinetics. The difference between these two is that Hill kinetics takes the cooperativity between proteins into account.8 The transcription and translation process is supported by multiple cooperative transcription factors, hence why Hill kinetics should be used for these processes. However, in this model, it is assumed that the only major transcription factor is dimerized STAT3 in the nucleus, and this transcription factor has no cooperativity with other transcription factors, resulting in a Hill coefficient of 1. The Hill equation with a Hill coefficient of 1, which is the case for transcription and translation, equals the Michaelis-Menten equation.8
All reaction constants were found in the literature and can be found in the parameters accordion. The dissociation rate constants (kd) and the extracellular degradation of IL-10 (KDe1) were determined with the dissociation constants, association rate constants, and half-life time. The dissociation rate constants were derived from the association constants (ka) and the dissociation constants (KD) with the formula kd = KD * ka. The degradation of extracellular IL-10 was determined from the half-life time of IL-10, which is 1 hour in vitro.9 IL-10 undergoes first-order kinetic elimination.10 So, the degradation rate constant (KDe1) can be determined with the formula KDe1 = ln(2) / t1/2.
Most of the initial concentrations are based on previously developed models by Yamada et al. and Moya et al.1,11 The species that are not present when a simulation starts, such as the complexes, mRNA, and IL-10, have a concentration of 0 nM. The ANCA concentration is based on the concentrations used in the laboratory and the concentration ANCAs present in ill patients. This is in the range of 10-40 ng/mL.12 The receptor concentration starts at 16.9 nM, but is changed during determining the optimal ratio between receptor and ligand.
Assumptions
The above-mentioned model is based on mass-action kinetics. For this, several assumptions are made. First of all, it is assumed that the cell is divided into two compartments (the cytoplasm and the nucleus) of equal volumes. Moreover, the volumes of the ER and Golgi complex, which are included in the model, are assumed to have equal volumes. Consequently, the concentration balance is maintained, allowing the transport from the ER via the Golgi complex to the extracellular environment to be simulated by mass-action kinetics. Furthermore, it is assumed that no ligand depletion in the extracellular environment occurs because of the considerably larger volume of the extracellular environment in comparison to the volume of the cell. As a result, the number of receptors is limited in comparison to the number of ligands.13 So, due to the large extracellular volume, ligand-receptor binding does not affect the concentration of the ligand.14 Therefore, the ordinary differential equation of the ligand concentration is set to zero, resulting in a constant ligand concentration over time.
| Abbreviation | Description |
|---|---|
| IL10er | IL-10 in Endoplasmic Reticulum |
| IL10ex | IL-10 in the extracellular |
| IL10g | IL-10 in Golgi complex |
| L | Ligand, ANCAs |
| mRNAc_IL10 | IL-10 mRNA in the cytoplasm |
| mRNAn_IL10 | IL-10 mRNA in the nucleus |
| PPN | Phosphatase in the nucleus |
| PPN_STAT3np | Phosphatase in nucleus bound to phosphorylated STAT3 |
| PPX | Phosphatase in the cytoplasm |
| PPX_STAT3cp | Phosphatase in the cytoplasm bound to phosphorylated STAT3 |
| RJ | Receptor with JAK bound |
| RJL | Receptor with JAK and ligand-bound |
| RJ2L | Dimerized receptor with JAK and ligand-bound |
| RJ2Lp | Dimerized and phosphorylated receptor with JAK and ligand-bound |
| RJ2Lp_STAT3 | Unphosphorylated STAT3 bound to the dimerized receptor |
| RJ2Lp_STAT3p | Phosphorylated STAT3 bound to the dimerized receptor |
| STAT3c | Unphosphorylated STAT3 in the cytoplasm |
| STAT3cp | Phosphorylated STAT3 |
| STAT3cpd | Dimerized phosphorylated STAT3 in the cytoplasm |
| STAT3n | Unphosphorylated STAT3 in the nucleus |
| STAT3np | Phosphorylated STAT3 in the nucleus |
| STAT3npd | Dimerized phosphorylated STAT3 in the nucleus |
| Name | Abbreviation | Initial value | Unit | Source |
|---|---|---|---|---|
| y1 | STAT3c | 1000 | nM | (Yamada et al., 2003) |
| y2 | RJ2Lp_STAT3c | 0 | nM | (Moya et al., 2011) |
| y3 | STAT3cp | 0 | nM | (Moya et al., 2011) |
| y4 | RJ2Lp_STAT3cp | 0 | nM | (Moya et al., 2011) |
| y5 | PPX | 50 | nM | (Yamada et al., 2003) |
| y6 | PPX_STAT3cp | 0 | nM | (Moya et al., 2011) |
| y7 | STAT3cpd | 0 | nM | (Moya et al., 2011) |
| y8 | STAT3npd | 0 | nM | (Moya et al., 2011) |
| y9 | STAT3np | 0 | nM | (Moya et al., 2011) |
| y10 | PPN | 60 | nM | (Yamada et al., 2003) |
| y11 | PPN_STAT3np | 0 | nM | (Moya et al., 2011) |
| y12 | STAT3n | 0 | nM | (Moya et al., 2011) |
| y13 | mRNAn_IL10 | 0 | nM | (Moya et al., 2011) |
| y14 | mRNAc_IL10 | 0 | nM | (Moya et al., 2011) |
| y15 | IL10er | 0 | nM | (Moya et al., 2011) |
| y16 | IL10g | 0 | nM | (Moya et al., 2011) |
| y17 | IL10ex | 0 | nM | (Moya et al., 2011) |
| y18 | L | 24 | nM | (Majid et al., 2021) |
| y19 | RJ | 16.9 | nM | (Reeh et al., 2019) |
| y20 | RJL | 0 | nM | (Moya et al., 2011) |
| y21 | RJ2L | 0 | nM | (Moya et al., 2011) |
| y22 | RJ2Lp | 0 | nM | (Moya et al., 2011) |
| Parameter | Original value | Lower boundary | Upper boundary | Unit |
|---|---|---|---|---|
| KD1 | 0.741 | 0.001 | 1000 | nM |
| KB7 | 0.0126 | 0.000126 | 1.26 | nM-1s-1 |
| KBu11 | 0.0093366 | 0.000000126 | 1260 | s-1 |
| Abbreviation | Description | Value | Unit | Source |
|---|---|---|---|---|
| KB1 | Rate of binding (ka) of STAT with receptor | 0.008 | nM-1s-1 | (Yamada et al., 2003) |
| KB2 | Rate of binding (ka) of STATp with receptor | 0.005 | nM-1s-1 | (Yamada et al., 2003) |
| KB3 | Rate of dimerization (ka) of STATp in the cytoplasm | 0.020 | nM-1s-1 | (Moya et al., 2011) |
| KB4 | Rate of binding (ka) of STATp with PPX | 0.001 | nM-1s-1 | (Yamada et al., 2003) |
| KB5 | Rate of dimerization (ka) of STATp in the nucleus (=KB3) | 0.020 | nM-1s-1 | (Moya et al., 2011) |
| KB6 | Rate of binding (ka) of PPN and STATp | 0.001 | nM-1s-1 | (Yamada et al., 2003) |
| KB7 | Rate of binding (ka) of the ligand with the receptor | 0.0126 | nM-1s-1 | (Granel et al., 2020) |
| KB8 | Rate of dimerization (ka) of RJL and RJ | 0.000001-0.0003 | nM-1s-1 | (Pang & Zhou, 2012) |
| KBu1 | Rate of dissociation (kd) of STAT and receptor | 0.8 | s-1 | (Yamada et al., 2003) |
| KBu2 | Rate of dissociation (kd) of STAT and receptor with phosphorylation | 0.4 | s-1 | (Yamada et al., 2003) |
| KBu4 | Rate of dissociation (kd) of STATp and receptor | 0.5 | s-1 | (Yamada et al., 2003) |
| KBu5 | Rate of dedimerization (kd) of STATp in the cytoplasm | 0.1 | s-1 | (Moya et al., 2011) |
| KBu6 | Rate of dissociation (kd) of STAT and PPX without dephosphorylation | 0.2 | s-1 | (Yamada et al., 2003) |
| KBu7 | Dephosphorylation of STAT by dissociation of PPX | 0.003 | s-1 | (Yamada et al., 2003) |
| KBu8 | Rate of dedimerization (kd) of STATp in the nucleus (=KBu5) | 0.1 | s-1 | (Moya et al., 2011) |
| KBu9 | Rate of dissociation (kd) of STATp and PPN without dephosphorylation | 0.2 | s-1 | (Yamada et al., 2003) |
| KBu10 | Rate of dissociation (kd) of STATp and PPN with dephosphorylation | 0.005 | s-1 | (Yamada et al., 2003) |
| KBu11 | rate of dissociation (kd) L and RJ, calculated with KD1 and KB7 | Calculate | s-1 | Calculated |
| KBu12 | Rate of dissociation (kd) of the two receptors, calculated with KD2 and KB8 | Calculate | s-1 | Calculated |
| KD1 | Range of values for dissociation constant of ANCAs with PR3 | 0.001 - 1000 | nM | (Granel et al., 2020) |
| KD2 | Dissociation constant of dimerization of the receptors | 1000 – 100000 | nM | (Philo et al., 1996) |
| KDn1 | Degradation rate in the nucleus | 0 | s-1 | - |
| KDc1 | Degradation rate in the cytoplasm | 0.0104 | s-1 | (Maiti et al., 2015) |
| KDe1 | Degradation rate in the extracellular | 0.0001925 | s-1 | Calculated |
| KM1 | Maximal velocity of transcription of IL-10 | 0.01 | nMs-1 | (Yamada et al., 2003) |
| KM2 | Michaelis-Menten constant for transcription of IL-10 | 400 | nM | (Yamada et al., 2003) |
| KM3 | Maximal velocity of translation of IL-10 | 0.006 | s-1 | (Marshall & Noireaux, 2019) |
| KM4 | Michaelis-Menten constant for translation of IL-10 | 10 | nM | (Marshall & Noireaux, 2019) |
| KMH1 | The Hill coefficient in the transcription of IL-10 | 1 | - | - |
| KP1 | Phosphorylation rate of the receptor | 0.005 | s-1 | (Yamada et al., 2003) |
| KTcSTAT3 | The transportation rate of STAT3 from the nucleus to the cytoplasm | 0.05 | s-1 | (Yamada et al., 2003) |
| KTnSTAT3cpd | The transportation rate of STAT3pd from the cytoplasm to the nucleus | 0.005 | s-1 | (Yamada et al., 2003) |
| KTcmRNA_IL10 | The transportation rate of IL-10 mRNA from the nucleus to the cytoplasm | 0.001 | s-1 | (Yamada et al., 2003) |
| KTerIL10 | The transportation rate of IL-10 from the Endoplasmic Reticulum to the Golgi complex | 0.0005 | s-1 | (Hirschberg et al., 1998) |
| KTgIL10 | The transportation rate of IL-10 from the Golgi complex to the extracellular | 0.0005 | s-1 | (Hirschberg et al., 1998) |
| Rtot | Concentration of ribosomes within a cell | 1100 | nM | (Marshall & Noireaux, 2019) |
| t12_iv | Half-life time of IL-10 in the extracellular | 3600 | s-1 | (Le et al., 1997; Milligan et al., 2005) |
Ordinary Differential Equations
Chemical equations
Multiple Parameter Sensitivity Analysis (MPSA)
Due to the large size and complexity of the model, relating parameters to their direct effects was challenging and could be easily compensated by other parameters. A multiple parameter sensitivity analysis (MPSA), as proposed by Zi et al., was performed on three parameters.15 These parameters were selected based on their adjustability in the lab. This way, the results of the MPSA can directly be used in the laboratory.
The sensitivity of the parameters related to one feature was determined, the IL-10 concentration at time = 36 hours. The feature is based on the read-out time in the laboratory, after around 36 hours. The parameters were sampled from a biologically feasible order of magnitude, and the number of parameter subset samples was investigated. The elaborate method that was used for this analysis can be found in the accordion below.
The parameters that are included in the MPSA are selected based on their ability to be adjustable in the laboratory. We suspect that changing the downstream parameters will have a greater influence on the IL-10 concentration over time than changing the upstream parameters. Whereas a signal pathway undergoes amplification, a change in the rate of the biochemical reactions at the end of the pathway will have a greater influence on the IL-10 production than a reaction at the beginning. From the results of the MPSA, the optimal values for the parameters are chosen.
Ten thousand parameter sets (N) were chosen uniformly using Latin Hypercube sampling (LHS), a statistical method generating near-random samples, where the set of random numbers is representative of the real variability.15 Parameters were logarithmic sampled through a range between a predefined upper and lower boundary. These boundaries were selected based on a biological relevant order of magnitude and can be found in appendix conditions.
The sensitivity (S) per feature (f) and set (I) was described based on the difference between a reference state y¬_ref, and the model output y_sim. For y_sim a preselected parameter set Θ(i),i=1,…, N, and for y_ref the Il-10 reference set was used, see Equation 1.
S(f,i) = (yref(f) – ysim(Θ(i,f))2 Equation 1
Two classes were created based on a threshold value. The threshold value was determined, based on the feature mean value. By comparing the feature values of the sets with the threshold a difference can be seen between parameter sets with a large or small differences in reference to the mean see Equation 2.
After classifying N parameter sets, the cumulative frequency function of the two classes was determined per parameter (i).
S(f,Θ(i)) ≥ Treshold Θ(i)∈Set Shigh(f) Equation 2a
S(f,Θ(i)) < Treshold Θ(i)∈Set Slow(f) Equation 2b
The Kolmogorov-Smirnov score (D) was used to quantify the sensitivity of each feature f. Here D is the maximal distance between the cumulative frequency functions, and an increasing D represents an increase in sensitivity, see Equation 3.
D = sup(Θ) | Sa(Θ) – Su(Θ) | Equation 3
Ten dummy parameters were added, and different sample sizes N were tried to determine the minimal sample size for the MPSA. Dummy parameters do not influence the output of the model but have a non-zero sensitivity index. This error is due to the numerical approximation and, therefore, an often-used measure to determine the MPSA sample size.
- Yamada S, Shiono S, Joo A, Yoshimura A. Control mechanism of JAK/STAT signal transduction pathway. FEBS Lett. 2003;534(1-3). doi:10.1016/S0014-5793(02)03842-5
- Singh A, Jayaraman A, Hahn J. Modeling regulatory mechanisms in IL-6 signal transduction in hepatocytes. Biotechnol Bioeng. 2006;95(5). doi:10.1002/bit.21026
- Heinrich PC, Behrmann I, Haan S, Hermanns HM, Müller-Newen G, Schaper F. Principles of interleukin (IL)-6-type cytokine signalling and its regulation. Biochemical Journal. 2003;374(1). doi:10.1042/BJ20030407
- Lokau J, Schoeder V, Haybaeck J, Garbers C. Jak-stat signaling induced by interleukin-6 family cytokines in hepatocellular carcinoma. Cancers (Basel). 2019;11(11). doi:10.3390/cancers11111704
- Stahl N, Boulton TG, Farruggella T, et al. Association and activation of Jak-Tyk kinases by CNTF-LIF-OSM-IL-6 β receptor components. Science (1979). 1994;263(5143). doi:10.1126/science.8272873
- Kojima R, Aubel D, Fussenegger M. Building sophisticated sensors of extracellular cues that enable mammalian cells to work as “doctors” in the body. Cellular and Molecular Life Sciences. 2020;77(18). doi:10.1007/s00018-020-03486-y
- Kim M, Morales LD, Jang IS, Cho YY, Kim DJ. Protein tyrosine phosphatases as potential regulators of STAT3 signaling. Int J Mol Sci. 2018;19(9). doi:10.3390/ijms19092708
- Singh V. Introduction to synthetic biology. In: Advances in Synthetic Biology. ; 2020. doi:10.1007/978-981-15-0081-7_1
- Le T, Leung L, Carroll WL, Schibler KR. Regulation of interleukin-10 gene expression: Possible mechanisms accounting for its upregulation and for maturational differences in its expression by blood mononuclear cells. Blood. 1997;89(11). doi:10.1182/blood.v89.11.4112
- Milligan ED, Langer SJ, Sloane EM, et al. Controlling pathological pain by adenovirally driven spinal production of the anti-inflammatory cytokine, interleukin-10. European Journal of Neuroscience. 2005;21(8). doi:10.1111/j.1460-9568.2005.04057.x
- Moya C, Huang Z, Cheng P, Jayaraman A, Hahn J. Investigation of IL-6 and IL-10 signalling via mathematical modelling. IET Syst Biol. 2011;5(1). doi:10.1049/iet-syb.2009.0060
- Reeh H, Rudolph N, Billing U, et al. Response to IL-6 trans- A nd IL-6 classic signalling is determined by the ratio of the IL-6 receptor α to gp130 expression: Fusing experimental insights and dynamic modelling. Cell Communication and Signaling. 2019;17(1). doi:10.1186/s12964-019-0356-0
- Finlay DB, Duffull SB, Glass M. 100 years of modelling ligand–receptor binding and response: A focus on GPCRs. Br J Pharmacol. 2020;177(7). doi:10.1111/bph.14988
- Webster TF, Schlezinger JJ. Generalized concentration addition for ligands that bind to homodimers. Math Biosci. 2019;316. doi:10.1016/j.mbs.2019.108214
- Zi Z. Sensitivity analysis approaches applied to systems biology models. IET Syst Biol. 2011;5(6):336-346. doi:10.1049/iet-syb.2011.0015
- Majid, T., Al-Obaidy, M., Ali Mutar, H., & Salih, T. A. (2021). Studying of Some Immunological Parameters in Patients with Inflammatory Bowel Disease ( IBD ) in Al-Anbar Province. Journal of University of Anbar for Pure Science (JUAPS) Open Access, 2021(1), 1–7. https://doi.org/10.37652/juaps.2022.172396
- Granel, J., Lemoine, R., Morello, E., Gallais, Y., Mariot, J., Drapeau, M., Musnier, A., Poupon, A., Pugnière, M., Seren, S., Nouar, D., Gouilleux-Gruart, V., Watier, H., Korkmaz, B., & Hoarau, C. (2020). 4C3 Human Monoclonal Antibody: A Proof of Concept for Non-pathogenic Proteinase 3 Anti-neutrophil Cytoplasmic Antibodies in Granulomatosis With Polyangiitis. Frontiers in Immunology, 11. https://doi.org/10.3389/fimmu.2020.573040
- Pang, X., & Zhou, H. X. (2012). A common model for cytokine receptor activation: Combined scissor-like rotation and self-rotation of receptor dimer induced by class I cytokine. PLoS Computational Biology, 8(3). https://doi.org/10.1371/journal.pcbi.1002427
- Philo, J. S., Aoki, K. H., Arakawa, T., Narhi, L. O., & Wen, J. (1996). Dimerization of the extracellular domain of the erythropoietin (EPO) receptor by EPO: One high-affinity and one low-affinity interaction. Biochemistry, 35(5). https://doi.org/10.1021/bi9524272
- Maiti, S., Dai, W., Alaniz, R. C., Hahn, J., & Jayaraman, A. (2015). Mathematical modeling of pro- and anti-inflammatory signaling in macrophages. Processes, 3(1). https://doi.org/10.3390/pr3010001
- Marshall, R., & Noireaux, V. (2019). Quantitative modeling of transcription and translation of an all-E. coli cell-free system. Scientific Reports, 9(1). https://doi.org/10.1038/s41598-019-48468-8
- Hirschberg, K., Miller, C. M., Ellenberg, J., Presley, J. F., Siggia, E. D., Phair, R. D., & Lippincott-Schwartz, J. (1998). Kinetic analysis of secretory protein traffic and characterization of Golgi to plasma membrane transport intermediates in living cells. Journal of Cell Biology, 143(6). https://doi.org/10.1083/jcb.143.6.1485
Model validation
MPSA
Parameter variations
Conclusion
Improvements
Outlook
References
Model validation
Internal correctness
Several tests were performed to determine whether the model satisfies biological feasible behavior. One of these requirements is that a species' total concentration, like STAT3, must remain stable over time. The receptor, PPX, and PPN concentrations should behave similarly. Another test is related to the amount of mRNA present, an increase of mRNA is seen when there is an absence of phosphatases or if there is lower STAT3-phosphatase binding. More mRNA ought to be produced if dimerized STAT3 in the nucleus is kept from disintegrating. Last but not least, the translation was examined to see if extracellular IL-10 was generated in the presence of mRNA. All the above-mentioned tests were successfully accomplished by the model, proving that it displays biological feasible behavior.
The system behavior
The model's behavior was evaluated using laboratory initial concentrations. The maximum anti-neutrophil cytoplasmic antibody (ANCA) concentration applied in the laboratory is 1000 pg/mL, which is equivalent to 6.6667 nM.
First, the receptor conversions are considered, as visualized in Figure 1. Continuous activation of the pathway by ANCAs resulted in an immediate drop in the GEMS receptor concentration. The rapid binding of the ligand to the receptor causes this drop. This binding converts the receptor into a receptor bound to the ligand, whose concentration immediately rises. Following the formation of a receptor monomer bound to the ligand, another receptor monomer can bind, resulting in the formation of a receptor homodimer. This reduces the concentration of monomer receptors and receptors bound to ligands, while increasing the concentration of dimerized receptors. However, because dimerized receptor phosphorylation is rapid and irreversible, no increase in dimerized receptors was observed. Nevertheless, there is an increase in dimerized and phosphorylated receptors. At first, there are many STAT3 proteins in the cytoplasm, which causes a rise in the concentration of STAT3 bound to the receptor and a decrease in the concentration of STAT3. When there is less STAT3 protein in the cytoplasm, the concentration of receptors bound to STAT3 decreases to an equilibrium.
STAT3 was phosphorylated as a result of its binding to the dimerized receptor. The concentration of the formed STAT3cp increased rapidly and reaches a steady state after 3 hours (Figure 2). Braun et al. (2013) observed a rapid increase in the concentration of the phosphorylated STAT3 in the cytoplasm, followed by a decrease until a steady state was reached. However, in their study, a steady state was obtained after 20 minutes.1
(B) Zoomed in on reaching a steady state of phosphorylated and dimerized STAT3 in the cytoplasm.
After the dimerization of two phosphorylated STAT3 monomers, the dimerized STAT3 translocated to the nucleus. The formation of dimerized STAT3 was faster than the dimer's breakdown and transport to the nucleus. Figure 3 shows an increase in dimerized STAT3. This rise is also due to an increase in the concentration of phosphorylated STAT3. Finally, the concentration decreases until a steady state is reached. After an increase in the concentration of dimerized STAT3 in the cytoplasm, the concentration of dimerized STAT3 in the nucleus increases and reaches a steady state after 3 hours. According to Yamada et al., the model shows an accumulation of the STAT3 dimer in the nucleus.2
Only a small amount of mRNA was present in the system due to cytoplasmic mRNA degradation (Figure 4). The Endoplasmic Reticulum (ER) produced the first interleukin 10 proteins after half an hour. When compared to the proteins from the ER, IL-10 proteins from the Golgi complex were produced with a delay. This makes sense given the time it takes to get ER to the Golgi complex. The rate constants for export from the ER and the Golgi complex are the same, so it makes sense that the steady-state concentration of IL-10 protein from the ER and Golgi complex have the same concentrations. Finally, the extracellular fluid IL-10 protein concentration is increased. After some time, this concentration also reaches a steady state because the rate of production and degradation of extracellular IL-10 becomes the same.
Multiple Parameter Sensitivity Analysis (MPSA)
The multiple parameter sensitivity analysis was used to determine which parameter is the most sensitive to the expression of IL-10. The analysis included KD1, KB7, and KBu11, all of which have the ability to be altered in the laboratory. KD1 is the dissociation constant between the GEMS receptor and the ANCAs, KB7 is the binding rate and KBu11 is the unbinding rate between the GEMS receptor and the ANCAs. The other parameters present in the ODE model are not easily adjustable in the laboratory and are, therefore, not considered. The dimerization between receptors we considered because a point mutation could be implanted and in this way the rate of binding and unbinding could be changed. It is only uncertain what kind of point mutations has what kind of consequences, and that is why they are left out of the MPSA.
The MPSA is performed with a sample size of 10000 parameter sets and ten dummies. The ten dummies were included to determine the sufficient sample size. The sample size should be high enough to get a minimal differentiating number for the sensitivity of the dummies, if this is not the case the sample size should be increased. Since the sensitivities resulting from the MPSA are then dependent on the sample values selected.
The sensitivities are determined for sample sizes logarithmically chosen, ranging from 500 to 10000 samples. Figure 5 shows that the sensitivity decreases with an increasing sample size increases, this was expected. Additionally, Figure 5 demonstrates that a sample size of 10000 is sufficient for the study because, at this level, the sensitivity does not decline anymore.
We expected that parameters downward in the pathway would have a higher sensitivity towards the interleukin 10 (IL-10) concentration than the upward parameters because of the amplification effect within a pathway. For instance, when a reaction downwards works slower than the whole pathway and all its amplifications will be inhibited until this point. The model is looked at, at the IL-10 concentration after around 36 hours, this is the same time as the read-out in the laboratory during an experiment.
In Figure 6 you can see the cumulative probability distributions of the unaccepted (green) and accepted parameter (yellow) sets (A, B, C) and the calculated sensitivities (D). The sensitivity is based on the maximal difference between the two cumulative distributions, see the methods section for the elaborated methods. A large difference between the two distributions visualizes that the corresponding parameter has a large influence on the classification of these two sets.
KD1 and KB7 are less sensitive than KBu11 by looking at Figure 6D. The parameter ranges could explain this difference in sensitivity. The parameter ranges are based on biological feasible ranges for all three parameters. KBu11 can be changed to much larger numbers than both KD1 and KB7. As a result, this parameter can influence the IL-10 concentration more.
Parameter variations
Receptor-ligand ratio
In patients with antineutrophil cytoplasmic antibody (ANCA)-associated vasculitis (AAV) is the assumed concentration of ANCAs present in the body around 35.78 ng/mL, this is equivalent to 0.24 nM.3 This initial ligand concentration is not changed because this is the amount present in patients. As a result, when more or less interleukin 10 (IL-10) should be produced by the synthetic cell, the concentration of the GEMS receptor on the surface of the cell should be changed. From the results of the general model, it was concluded that the optimal ratio of receptor-ligand would be around 100:1. When this ratio is applied to this model, the ligand's initial concentration is 0.24 nM and the receptor's initial concentration is 24 nM. Based on these initial concentrations, the receptor concentration was varied over 24*10-3 nM to 24*103 nM. In comparison to the ligand concentration, this ranges to a deficit or excess of the receptor.
Figure 7 reveals two intriguing values. (1) The minimal concentration of receptor required to form the first dimerized receptors and activate the pathway. (2) Starting from what concentration is the maximum RJ2Lp reached and should thus be the highest receptor concentration because any higher concentration is useless.
When the RJ2Lp begins to rise, the pathway is activated, and the receptor concentration is at its lowest. Lower receptor concentrations are ineffective because no RJ2Lp is formed, and the pathway is not activated. As shown in the graph (Figure 7), the receptor-ligand ratio is around 5 when the RJ2Lp concentration starts to rise. This means that at least 5 receptors must be present for each ligand to result in a dimerized receptor and thus activation of the pathway.
When adding more receptors does not affect the RJ2Lp concentration, the maximum receptor concentration is reached. The flattening of the curve at the high RJ2Lp concentrations indicates this. When a ratio of 500 is reached, the curve flattens and the maximum RJ2Lp concentration is reached. A ratio of 500 indicates that 500 receptors must be present for each ligand to achieve maximum RJ2Lp concentration; a higher receptor-ligand ratio is not required. The following simulations are conducted with a receptor-ligand ratio of 500.
Changing dissociation constant
The dissociation constant (KD) between the GEMS receptor and the ligand is the most easily adjusted parameter in the laboratory. The ligand is an ANCA, which is an antibody. Antibodies have a wide range of dissociation constants spanning from 10-6 to 10-12 M. Since the KD is antibody and GEMS receptor dependent, when a change in KD is required, another antibody can be ordered, or a change in the GEMS receptor can be made through a point mutation or by creating a different extracellular affinity domain.
Figure 8 depicts the IL-10 concentration over the range of KD values at the time point 36 hours. The maximum IL-10 concentration is reached after 36 hours for all KD values. It appears that the KD has little effect on the IL-10 concentration.
The impact of the KD1 can be visualized over time. For all KD1 values, an equilibrium of IL-10 concentration is reached after 36 hours. The maximum IL-10 concentrations are probably not reached at the same time, but this is not depicted in Figure 8.
Figure 9 illustrates the influence of the KD1 over time. Early on, KD1 values with a value closer to the picomolar range have a higher rise in IL-10 concentration. While KD1 values in the micromolar range have a slow increase in IL-10 concentration.
Conclusion
All the tests performed to determine whether the model satisfied biologically feasible behavior were successfully accomplished. The concentrations of the receptor, PPX, PPN, and STAT3 remained stable over time. The mRNA reacted as predicted with the phosphatases, and when mRNA is present, extracellular IL-10 is produced. This demonstrates that the model exhibits biologically feasible behavior. The model achieves equilibrium, which is desirable because almost all biological processes achieve a steady state for their concentrations.
For the MPSA, we expected that KD1 or KBu11 would have the biggest influence on IL-10 secretion because they influence the binding time. This was the case for KBu11, with a sensitivity of 0.2749, but it was not the case for KD1, with a sensitivity of 0.02303. KBu11 is the unbinding rate of the ANCAs from the receptor, by slowing down this rate the IL-10 production can be increased and the other way around. If the unbinding rate is too fast, the pathway cannot be activated and no IL-10 is produced. The laboratory should focus on the duration of the binding between ANCAs and receptors if the secretion of IL-10 needs to be influenced or optimized.
A maximum receptor-ligand ratio of 500 was deduced from the graph in Figure 7. This is significantly higher than expected from the previous model's maximum ratio of 100. Increasing the receptor concentration to an even higher receptor-ligand ratio does not affect the IL-10 concentration and is therefore ineffective. The minimum effective receptor-ligand ratio had to be at least 5. This is also higher than expected because a dimerized receptor consists of only two receptor monomers for each ligand, so we expected the minimum ratio to be around 2. A ratio of 5 is not remarkably higher, so this could be possible.
The MPSA results showed that varying the KD1 value at 36 hours had no effect on IL-10 concentration, which is also supported by Figure 8. The maximum IL-10 concentration is reached for all these KD1 because it is an equilibrium concentration. When the IL-10 concentration was examined at different time points for the range of KD1 values, it was discovered that a picomolar range causes a faster increase in IL-10 concentration and a micromolar range causes a slower increase in IL-10. Accordingly, when a high concentration of IL-10 is required in a very short period of time, highly specific antibodies for the affinity domain are required. To summarize, KD1 affects the IL-10 concentration, but the effect is time-dependent.
Improvement
In comparison to the literature, the time it takes to reach steady states is quite slow.1,2 This could probably be solved by lowering the dissociation constant of receptor dimerization (KD2), which would result in faster signaling. Therefore, re-estimating the values of KD2 would be a suggestion for future research. At the moment, the dimerization dissociation constant is based EpoR dimerization, whereas only the affinity domains dimerize. Selecting a different affinity domain to dimerize with, could result in faster gene expression.
Outlook
The model could be extended by adding a feedback loop between the IL-10 production and the ANCA concentration. This has not yet been done because the relationship between IL-10 and ANCAs is unknown. Nonetheless, such a feedback loop could provide insight into the progression of the disease as well as the effect of our synthetic cell-based treatment. Adding the feedback loop does not make sense for our treatment because the relationship between the treatment and disease is unknown, but the model can be used by others, so if the feedback loop is known for their treatment, this would be a great way to extend the model.
- Braun DA, Fribourg M, Sealfon SC. Cytokine response is determined by duration of receptor and signal transducers and activators of transcription 3 (STAT3) activation. Journal of Biological Chemistry. 2013;288(5). doi:10.1074/jbc.M112.386573
- Yamada S, Shiono S, Joo A, Yoshimura A. Control mechanism of JAK/STAT signal transduction pathway. FEBS Lett. 2003;534(1-3). doi:10.1016/S0014-5793(02)03842-5
- Majid T, Al-Obaidy M, Ali Mutar H, Salih TA. Studying of Some Immunological Parameters in Patients with Inflammatory Bowel Disease ( IBD ) in Al-Anbar Province. Journal of University of Anbar for Pure Science (JUAPS) Open Access. 2021;2021(1):1-7. doi:10.37652/juaps.2022.172396