Modularisation of published and novel models toward a complex KIR2DL4 pathway in pbNK cell

KIR2DL4 is an interesting receptor expressed on the peripheral blood natural killer (pbNK) cell as it can be either activating or inhibitory depending on the amino acid residues in the domain. This model uses mathematical modelling to investigate the downstream effects of natural killer cells’ activation (KIR2DL4) receptor after stimulation by key ligand (HLA-G) on pbNK cells. Development of this large pathway is based on a comprehensive qualitative description of pbNKs’ intracellular signalling pathways leading to chemokine and cytotoxin secretion, obtained from the KEGG database (https://www.genome.jp/pathway/hsa04650). From this qualitative description we built a quantitative model for the pathway, reusing existing curated models where possible and implementing new models as needed. This model employs a composite approach for generating modular models. The approach allows for the construction of large-scale complex model by combining component of sub-models that can be modified individually. This large pathway consists of two published sub-models; the Ca2+ model and the NFAT model, and a newly built FCεRIγ sub-model. The full pathway was fitted to published dataset and fitted well to one of two secreted cytokines. The model can be used to predict the production of IFNγ and TNFα cytokines.• Development of pathway and mathematical model• Reusing existing curated models and implementing new models• Model optimization and analysis


Background
Natural killer (NK) cells being primarily associated with innate immunity. There is some evidence that they also contribute to adaptive immunity, especially in adapting the environment through priming, education and memory, and are able to respond to multiple infections from a particular antigen [1,2,3]. They comprise approximately 15% of all circulating lymphocytes [4,5,6] including 2-18% of lymphocytes in human peripheral blood. NK cells have the capability to lyse or destroy some abnormal cells -for example, tumour cells and virus-infected cells and are known as effector cells, which means they selectively bind to ligands and regulate target cells' biological activities [7,8]. NK cell function is determined by a balance of activation and inhibition signalling induced by trans-membrane receptors [9,1]. This process involves engagement of ligands with the receptors, as well as the action of pro-inflammatory cytokines such as IL-1, IL-2, IL-12, IL-15, IL-18, IL-21 and IFNα, β [7,10,11]. If activation signalling dominates, NK cells become activated. Upon contact with other cells recognised as undesirable, activated NK cells are able to mediate cell killing via two mechanisms, exocytosis of perforin/granzyme granules and signalling via TNF receptors [12].
The interaction between NK cells and target cells requires docking between NK cell receptors and target cell ligands. In this study, we developed a model to predict the production of IFNγ and TNFα cytokines induced by the binding of HLA-G to endocytosed KIR2DL4 receptor as shown in Figure 1 The role of PLCγ is to control the concentration of calcium in the cell [13,14] through IP3. PLCγ catalyses the conversion of IP3 from PIP2. IP3 causes the release of Ca 2+ from the intracellular endoplasmic reticulum. Ca 2+ binds to calcineurin (CaN), and the complex then dephosphorylates the nuclear factor in an activated T-cell (NFAT).
The migration of NFAT into the cell nucleus starts the transcription of cytokines TNFα and IFNγ [15,16].
The signalling pathway to be modelled, depicted in Figure 2 consist of many components and reactions. We apply modular design principles into the construction of NK signalling model. This model incorporates two published models including the Ca 2+ oscillations of Dupont and Erneux (1997) [17] and NFAT cycling of Cooling et al. (2009) [18]. A new sub-model was built to complete the signalling pathway, based on the availability of experimental data. One sub-model was developed in this section and termed FCRIγ and appear as Level 2 models in Figure 2. With all required sub-models developed, we were then able to bring together these models into a complete HLA-G-KIRDL4 pathway and labelled as Level 3 in Figure 2.

Results
Re-using and implementing models from Physiome Model Repository (level 1 models) Ca 2+ oscillations in Dupont and Erneux, 1997 This section outlines the model labelled as Ca 2+ sub-model in Figure 2. Dupont and Erneux simulated a model that captures the oscillation of Ca 2+ in a cell in response to upstream signalling ( Figure 3) [17]. The repository associated with this model can be found at https://github.com/Nurulizza/HLAG_to_cytokine/tree/master/dupont_1997.
Our pathway does not include the production of IP4 by phosphorylated IP3. To exclude this reaction, we zeroed the V3k reaction in Dupont's model which represents the activity of IP3 3-kinase in catalysing the IP4 production. In this case, when the maximal velocity of 3-kinase is zero, and the activity of 5-phosphatase very much exceeds the activity of 3-kinase, oscillations in IP3 disappear. The behaviour has also been observed in studies of other cell types [19,20,21]. The Ca 2+ and IP3 behaviour in the system without the activity of 3-kinase is illustrated in Figure 4. As expected, we retained the oscillation of Ca 2+ and the oscillation of IP3 disappeared without the changes of IP3 into IP4. Cooling et al., 2009 This section describes the model labelled as NFAT sub-model in Figure 2. Cooling et al. [18] adapted the experimental protocol of Tomida et al. [22], as described in previously to model NFAT cycling ( Figure 5) FCRIγ appears as Level 2 models in Figure 2. The repository associated with this model can be found at https://github.com/Nurulizza/HLAG_to_cytokine/blob/master/FCepsilonRI.cellml.

NFAT cycling in
The simulated Grb2 phosphorylation was a good fit with the experimental observations of Tsang et al. [23]. The model and observed experimental kinetics of phosphorylation of Grb2 are shown in Figure 7. The root mean square (RMS) error for the model is 0.0695 µM. The RMSE is 1.07% of the peak concentration. Visual inspection of the graph confirms this. The best fit to this data set was achieved with the parameters listed in Table 1.
Knowledge of which parameters the model is sensitive to helped to establish which parameters needed to be fixed when we tried to predict to another set of data (eg. the phosphorylation of Grb2 parameters (data not shown). The latter two types of parameter could therefore be varied over the whole defined range in fitting to subsequent data sets.
Parameter fitting was then performed by fixing the parameters that the previous fitted model was sensitive to: we fixed the value of kf2, kf4 kf6 and kf7 (see Table 2).
In this subsequent fitting, for constant forward rate parameters kf1 and reverse rate parameters kr1, kr4 and kr6 that showed distinct changes in model behaviour on one side of the minimum, we varied the parameters only over the range that model predictions for the Tsang data were not sensitive to these parameters. The list of kinetic rate constants that we are going to fit and the boundaries are given in Table   3.
The best fit was achieved with the parameters listed in Table 4. The model and the observed experimental [24] kinetics of phosphorylation of FC and Syk are shown in  Figure 8 these best fits are not perfect representations of the data, however, the available data are relatively sparse in time, and appear to contain more noise, so it may not be possible to obtain a clearly better fit to these data points. However, the model predictions shown in Figure 8 does not predict the experimentally observed drop in pSyk after 1000 s due to formation of pSyk-Grb2 complex. This may be because the concentration of Grb2 in the model was low (0.01 µM), which limited its impact on pSyk. The initial condition of Grb2 was fitted at 0.01 µM to get a good fit to the experimental data. The parameters for FCRIγ model fitted to data set from Tsang et al. [23] and Faeder et al. [24] are listed in Table 5. The sensitivity analysis of the model is provided in Additional File 2.
The sensitivity analysis showed that the model fitted to the data measured by In our model, the Syk is an OR-gate switch, meaning that it can reach full activation with one factor by ITAM binding. The ability of Syk to reach full activation with a single stimulus helps to define the ability of Syk to sustain its activity over time although after transient activation of ITAM to facilitate longer-term changes in cell signalling. As an example, Syk activity is required for more than 1 hour to induce activation of NFAT transcription [23].
Phosphorylation of Grb2 is dependent on the available concentration of Syk and In conclusion, a minimal FCRIγ model was developed and implemented in CellML. The results for phosphorylation of Grb2 for a data set from Tsang et al. [23] were presented to validate and test for phosphorylation of FCE and Syk in a data set from Faeder et al. [24]. A basic sensitivity analysis was performed to study the effect of the model constants' parameters on the model.

HLA-G cytokines model
This section introduces the model labelled HLA-G cytokines in Figure 2. The sensitivity analysis shows that IFNγ and TNFα productions are sensitive to both parameters, meaning that the changes of any parameter will affect the secretion of both cytokines. To elevate the secretion of TNFα at 8,000 s causes bad fit to IFNγ secretions. In this multi-objective optimisation, the model was only able to represent the IFNγ data accurately. When more data become available in the future, especially the knowledge of signalling pathway inside the nucleus, we hope the fitting can be improved. The model, however, captured the delay in cytokine production as expected. The best fit achieved with the parameters listed in Table 6.
The sensitivity analysis of the model is provided in Additional File 3. The sensitivity analysis of the model shows that the model is sensitive to parameters kf2, kf21, kf22, kf23 and kr21. The model is also sensitive to parameters kf4, kf5 and kf31 in component cytokines. The analysis showed that parameters kf21 depicts the binding activated calcineurin to NFATpc to form the complex NFATNc, kf21 and the nuclear import rate, kf22 are the most sensitive parameters that gave big difference to the model within the tested parameter space, that are between 10 0 and 10 2 and between 10 −2 and 10 2 . Parameters kf23 and kf2 sit on minimum within the the edge of regions within the tested parameters space. Parameters kf5 and kf31 are significant to either the production of IFNγ or TNFα.

Discussion
Here we aimed to model the regulation of activated KIR2DL4 receptor towards the production of IFNγ and TNFα cytokines. The activation of the endocytosed KIR2DL4 receptor was stimulated by the binding of soluble HLA-G. The full path-way was fitted to an NK cell dataset from Rajagopalan et al. [25]. As explained earlier, the use of the antibody KIR2DL4-specific involved small differences from natural HLA-G. In the model, we considered a simple case, which adequately catered for the biological understanding of the reactions. The model developed in this study has some limitations.
The sub-models are fitted to several types of data. We used B cell data [23] and RBL data [24] for the FCRIγ model and data for the pub-lished model by Dupount and Erneux [17] and Cooling et al. [18] were from Chinese hamster ovary cells and myocyte data, respectively.
In this study an ODE model has been employed. As the aim of the study is to capture a large signalling pathway, with significant interconnectivity, this provides the simple to generate and readily interpretable model of the system. In addition, information on spatial distribution of components of the pathway and data with which to parameterise a stochastic model is limited. However, in the future as data becomes available, it may be possible to derive more complete descriptions of the system using these techniques.
A J flux for each reaction is defined to simplify the equation. J flux comprises the concentration of species, and rate constants kf (forward reaction) and kr (reverse reaction). The model is encoded in CellML.
The signalling pathway to be modelled, depicted in Figure 2 consist of many components and reactions. To create a monolithic mathematical model which describes all these reactions would lead to a complicated model with many ODEs to solve. To simplify the problem, in this study, the development of the KIR2DL4 intracellular signalling model employs a composite approach for generating modular models.
The approach allows for the construction of large-scale complex model by combining component of sub-models that can be modified separately or individually [26]. The composite modular model brings together a series of sub-models. The sub-models can be derived from curated published models or newly developed models. (v) Sub-models may be independently calibrated or tested [26].
A major concern when modelling a biological system is that not all initial conditions and parameters for each activity in a cell are known. Unknown parameter values can be estimated using model calibration or fitting techniques to get values that fit the model behaviour to some experimentally observed behaviour. Model fitting can provide possible parameters that will help to produce reliable computational modelling results. In this study, the model parameter optimisation of CellML models was done using Python inside the OpenCOR tool. Model parameter optimisation to available experimental data helped us to simulate the model and capture experimental outputs.
In general, sensitivity analysis is the study of the changes of optimal solution when there are changes in the constant parameters [27]. A sensitivity analysis was performed to determine the sensitivity of rate constant parameters in the model. Each of these parameters determine changes in the model.

Conclusions
The The pathway was then fitted to pbNK cell cytokine secretion experimental data published by Rajagopalan et al. [25]. This approach was enabled through the availability of these models in standard formats in publicly accessible repositories. We apply modular design principles into the construction of signalling model.
The model was developed in CellML, which allows for modulation of the model. The model was developed in sub-models. The sub-models were developed so that they can be re-used in a new model, partially or as a whole. Each sub-models was fitted to available experimental data. An important contribution to the state of knowledge is the free availability of the code base, which makes it possible to independently simulate and to modify or extend the model for different applications.

Model construction
To identify appropriate signalling pathways to derive models at all levels and mod-els that had been previously published that covered components of the system we used online databases as resources for signalling pathway and experimental data. We then tested the model to ensure it worked as expected.
Methods for searching for experimental data Several methods were used to gain information for the experimental data in this study.
Various keywords for data searches were used, ranging from general to specific terms, Once a model has been developed, enough experimental data is crucial to ensure that the in silico simulation is reliable. The data sources can vary from journals to laboratory findings from experiments. Our model incorporated a wide variety of empirical observations ranging from human to non-human primates. Data that we needed included initial conditions of species involved in the pathway and kinetic rates of reactions. We examined literature data including integrative pathway diagrams, quantitative and qualitative details of possible molecular states, interaction and activities [28]. Appropriate parameter values were chosen from the literature where possible and others were optimised using parameter fitting.

Initial conditions
In this study, the initial conditions for substrates were derived from literature where possible. The initial conditions depend on the experiment conducted, or be fit if unknown. For experimental data that using the same unit, the value can be used directly. In some circumstance, experimental data was generated in different unit to the model, i.e. pg/ml, so the data need to be converted to molarity using a formula as follows

Molar concentration, = ci
where p is the density of constituent i, and M is the molar mass of constituent.
For substrate without any data, we used the relative initial condition that has been optimised for similar reaction or we made early assumption within the relevant range.
Once we optimised the model, we then decide whether we need to increase or decrease the value. The complexes in the system were fixed at 0 uM at time t =0.

Model calibration (parameter estimation)
Parameter estimation is a widely used method to estimate unknown parameters for mathematical models [29]. Most of the time, experimental data can't be used directly in a mathematical model due to many reasons, i.e. different experimental condition and measurement unit. Gutenkunst et al. highlighted that parameter estimation is complicated because of large measurement error and the model be-haviour also often insensitive to changes of an individual or combined parameter values [29,30].
Here we address the parameter estimation using least-squares optimisation in Python. The parameter estimation implements the Levenburg-Marquardt gradient method (greedy algorithm) to minimise an objective function [31]. The first step is to define the objective function to minimise. Gradient methods such as Levenburg-Marquardt tend to run into the nearest local minimum. The method is also sensitive to initialisation of parameters to be fit. The algorithm iteratively solves a trust-region sub-problems augmented by a special diagonal quadratic term and with trust-region shape determined by the distance from the bounds and the direction of the gradient.
This enhancements help to avoid making steps directly into bounds and efficiently explore the whole space of variables. Traditional gradient-based local optimization methods normally fail to get to a global minimum.
As the least-square method is prone to finding a local minimum rather than a global minimum [32,33], it is important to be reasonably close to the global minimum to get the best fit. To surmount this limitation, we must conduct a sweep of parameter space to find suitable search regions [34]. However, for many of our models, there are several unknown parameters, and so parameter space is too large to practically cover simply by varying parameters one by one. We therefore, aim to randomly sample parameter space with suitable coverage to be able to assure that we can find solutions close to a global minimum [35]. We use a sampling method to create a random sample of parameters within identified boundaries [36,37,38]. It allows users to determine the sampling size where the big sample size gives more chance for a better fit. It then runs the model lots of times using sampled parameter space and ranks the parameter set by the objective function. The parameter samples were generated using the Saltelli sample function stored in Python modules [36,37,38] as following where N is the scaling factor of samples to generate (the argument we supplied) and D is the number of model parameters. Parameter sweep generates samples within the specified parameter space.
We then generally loop over each sample input and evaluate the model. The Parameter estimation aims to find the possible parameters that minimise the difference between the model output and experimental data. We can measure how good is the fit by looking at the error value, where the lower error value means better fitting. We also can compare the graph between prediction and experimental data visually.
In this study, we used the default tolerance for termination by the change of the cost function and termination by the change of the independent variables. Default value set by the leastsq for both tolerances is 1e-8. The lower and upper bounds on independent variables are also set for each parameter.

Model analysis (sensitivity analysis)
Model uncertainty often has done as an analysis that is carried out after model construction and calibration have been completed [39]. One of the methodologies for uncertainty assessment is sensitivity analysis [39]. Sensitivity analysis is the study of how the inputs of a given model can change solution or output of the model, particularly of how the different level in the input of a model can be qualitatively or quantitatively apportioned to different outputs of a model [39,35,40]. One of the limitations of sensitivity analysis is the tendency to takes the model structure and system boundaries for granted [41].
There are many types of sensitivity analysis methods that can be used to access model uncertainty. For example, the quantitative variance-based methods, global sensitivity with regional properties [41], and the simplest class of the one factor at a time (OAT) screening techniques, which simply vary one factor at a time and measure the variation in the output [41]. Some form of various global sensitivity methods has been described in a few publications [35,42,43,44]. Ratto  be used to achieve a better understanding of parameter sensitivity [44,45].
The drawback of developing a big model in a small compartment is that the non-sensitive parameter from the small sub-model can be inherited to the bigger sub-model. Sensitivity analysis is one of the methods to access uncertainty in model parameter and can be used to identify parameter that need to be quantified experimentally.
To perform the analysis to NK cell model, each sub-model was analysed by looking at how sensitive is the model to each involved parameter. For small sub-models we use the one parameter at a time approach. That is, we change one parameter at a time and keep the other parameters fixed and investigate how that parameter impacts the objective function error. We initialised the model with the best fit values for each parameter, and ran it multiple times over the parameter range. We varied each parameter from its best fit value, over the parameter range we assumed it to take. The analysis is repeated for all the parameters. We vary parameters over the feasible parameter range to check how much they impact solutions.
Error defines the differences between points in model solution and experimental solution.
In this study, we used least-squares method which calculated the sum of the squared points of model from data. We plot a graph of the error (shown on the y-axis) versus parameter range (shown on the y-axis). Dips in the error correspond to local minima in the parameter space. If a model is sensitive to a parameter, it remained fixed in subsequent simulations.  The model also simulated the dephosphorylation of NFAT in the cytosol. From their model, they also found out that the dephosphorylated NFAT in the cytosol is buffered before declining. Nuclear NFAT was observed to rise steadily, as seen in the experimental data. The model was also able to replicate the effect of overexpression of calcineurin on the NFAT cycle. The model showed a normal NFAT cycle with overexpressed calcineurin at the same time that calcium is set at a constant level. In our cell signalling, the translocation of NFAT into the cell nucleus leads to secretion of cytotokines.
To simplify the model, as in the original publication [18], NFAT phosphorylation and dephosphorylation reactions are modelled as single steps, although in nature there is more than one phosphorylation site involved. The model also assumed a single step for NFAT translocation. The NFAT translocations into and from the nucleus do not have back reactions. This process is described using a futile cycle in the system.
With these assumptions, the phosphorylation and dephosphorylation of NFAT, together with translocations of NFAT into and from the nucleus, are described using mass action kinetics. The phosphorylation and dephosphorylation are described as where Nactive is the amount of activated calcineurin, determined by the total calcineurin, Ntot and fraction of activated calcineurin, ActN .
Thus, the translocations of NFAT into and from the nucleus are described as Parameter fitting was performed against the data of Tomida et al. [22], by Cooling [18] and the model was verified against the paper.

Implementing novel pathway models (level 2 models)
A new sub-model that was built to complete the whole signalling pathway, based on the availability of experimental data. One sub-model was developed in this section and termed FCRIγ and appear as Level 2 models in Figure 2.  [24], however, this model was focussed on the receptor itself and contained a greater level of complexity than was required in our model. In this sub-model we aimed to build a simple model of signalling of FC that still rep-resented biological understanding and fitted the experimental data ( Figure 6). FC exists in tetrameric form, which has an αchain, a β-chain and two disulfide-linked γ-chain [24]. The α subunit can be bound to ligand or unbound, a β-chain has four possible states of phosphorylation and two disulfide-linked γ-chain has six states of phosphorylations. In total, there are 300 possible states of receptor subunits for FC [24]. To avoid complexity, we assumed that FC is in an aggregated state and ready for binding with other reactants. The activated FC receptor is transphosphorylated by pLyn, available at the cell surface. This followed by the recruitment of Syk to the membrane surface which is subsequently transautophosphorylated by phosphorylated FCE. Inactive Grb2 then binds to phosphorylated Syk, and becomes phosphorylated.
Syk is a tyrosine kinase that is important in bridging receptor ligation and down-stream signalling such as Ca 2+ and MAPK. Once the cell receptor binds with the ligand, FCRIγ (ITAM receptor) is recruited and transphosphorylated by Lyn. The phosphorylated ITAM then recruits protein tyrosine kinase (Syk). Previous studies have shown that ITAM phosphorylation increases Syk activity and modulates Syk potency [23]. It must be noted that concentration above 10 µM ITAM can inhibit Syk activity [23]. Phosphorylated Syk then phosphorylates Grb2 (growth-factor-receptor-bound protein 2). Phosphorylation of Grb2 leads to the activation of PI3K downstream signalling. To match the available experimental data, we assumed that the Lyn kinase is rephosphorylated for recruitment of additional pLyn into the system to ensure that the system has enough supply of phosphorylated Lyn. The equation associated with this reaction is where Pi denoted a phosphate, k denotes a forward reaction and pLyn denotes a phosphorylated Lyn.
The flux associated with this reaction is where J denotes a flux, k denotes kinetic constant and Pi denotes a phosphate.
Experimental detail on Syk activity in B cells is presented in the study by Tsang et al. [23] which captures the molecular mechanism of Syk activation in vitro. The studies explore Syk activity under various conditions: activities of phosphorylated and nonphosphorylated Syk; activation of Syk by different receptors; and whether or not Syk is an OR-gate type of molecular switch. Syk activity is required for more than one hour to induce activation of downstream reactions, so the experiment was run for 3600 s. For that reason, the model simulation time was set for 3600 s. For the FCE model, we had 10 unknown parameters to fit to this data.
In one of the observations, Tsang looked at the activity of Syk when bound to FCRIγ.
The experimental data we fitted the model to were for a temporal Grb2 phosphorylation by dephosphorylated Syk with 1 µM FCRIγ. The Syk concentration used in the experiment was 0.005 µM. Experimental data points were digitised from Supplemental Figure 2 in Tsang et al. [23]. The observed experimental data points were from spectrophotometrical measurements of phosphorylation in a single representative experiment [23]. The binding of Syk to the receptor eliminated the lag phase in Grb2 phosphorylation.
The second data set was adapted from Faeder et al. [24]. In their simulation, Faeder et al.
looked at the pathway in rat basophilic leukemia (RBL) cells. The experimental environment was maintained at a temperature of 27 • C. The cell density and cell volume were assumed to be 1×10 6 cells/ml and 1.4×10 9 ml. The observed experimental data points were from densitometric measurements of phosphorylation in a single representative experiment [24,46]. Experimental data points were digitised from Figures 3(b) and 3(d) in Faeder et al. [24].
Parameter estimation was performed with each unknown kinetic rate constant allowed to vary in a feasible range that is between 1×10 −3 and 1×10 2 [47]. However, we expand the range when we see the optimisation error potentially decreasing with larger boundary. The settings for this work are described in Table 7. We estimated the concentration of Grb2 and pLyn at 6.47 µM and 6.5 µM. The phosphorylated (except Lyn) and complex reactants are assumed to be 0 at time 0 s. The number of sample used for the fitting is 500 which generate 12000 samples.
Tsang and Faeder used different experimental protocols so the initial conditions in our model representing these protocols differed. To fit the model to data from Faeder, some initial conditions were derived from the experimental protocol re-ported, and others were estimated during model fitting. We fitted concentration of FCRIγ is 0.0474 µM, the concentration of Syk is 0.025 µM, and the concentration of pLyn is 0.0474 µM. The list of initial conditions is listed in Table 8.

Modularisation of sub-models toward a complete HLA-G-KIRDL4 pathway
This section introduces the model labelled HLA-G cytokines in Figure 2. The respository associated with this model can be found at https://github.com/Nurulizza/HLAG_to_cytokine. Experimental studies by Rajagopalan et al. [25] observed the activation of KIR2DL4 by soluble HLA-G. The authors also suggest that the activation of uNK cells promotes a proinflammatory/proangiogenis response in uNK cells, which further promotes enlargement of blood vessels in early pregnancy [25].
Soluble HLA-G has been known to be a ligand for 2DL4. Soluble HLA-G originates from FCRIγ activation is also known to induce phosphoinositide 3-kinase (PI3K) [48,49].
Our model replicates the reaction of NK cells stimulated by KIR2DL4-specific antibody.
However, stimulation by the natural ligand, HLA-G, showed a quantitatively similar cytokine response. For most of the secreted cytokines, natural HLA-G induced at least 50% of the amount induced by antiKIR2DL4 mAb [25]. In the experiment, the secretion of TNFα was detected within the first 2 hours and up-regulated more than 2-fold. The secretion of IFNγ was upregulated 1.5-fold after 8 hours of receptor activation.
This large pathway is integrated by three level 1 and level 2 sub-models. Two of the submodels are published sub-models: the Ca 2+ model (labelled Ca 2+ in Figure 2) [17] and the NFAT model (labelled NFAT in Figure 2) [18]. This pathway also includes the FCRIγ model The fluxes associated with these reactions are Cytokine produces by a cell will be used for cell regulation. Although the model only captures the production of cytokine inside the cell, we added cytokines released from the cell that is applicable to all cytokines. We, therefore, added reactions for cytokine release for both TNFα and IFNγ.
The reactions associated with the reactions are Parameter fitting was performed using experimental data by Rajagopalan et al. [25].
Each unknown kinetic rate constant varies differently and the boundaries used for the fitting process are listed in Table 9. At this point, we fixed the initial conditions and kinetic rate constants estimated in individual published sub-models. However, we fitted parameters that were not sensitive to FC model.
We applied initial conditions as observed by Rajagopalan et al. [25] for HLA-G and KIR2Dl4. The initial conditions for HLA-G and KIR2Dl4 were set at 0.098 µM and 0.1052 µM, respectively. The initial conditions for PI3K was set at 0.01 µM. As shown by the reaction equation above, we assumed that the production of IFNγ and TNFα secretion were parallel to NFAT de-phosphorylation in the nucleus. We estimated the initial condition for plc (non-activated plc) was 1.3 µM, consisent with the initial condition for activated plc. The list of initial conditions adopted from the experiment for the HLA-G activation model, NFAT cytokines model and PI3K activation are shown in Table 10. In this model, we fitted 13 kinetic parameters including 5 parameters from the FC model.              Table 9 Parameters to fit to the data in Rajagopalan et al. [25] and the boundaries condition.

Parameter
Lower bound Upper bound   Figure 1 Schematic representation of TNFα and IFNγ secretions induced by HLA-G signalling pathway.

Figure 2
A phylogenetic tree showing the evolutionary relationships among existing and new sub-models of the pathway. Sub-models highlighted in green boxes mean experimental data exists.

Figure 4
Simulation of Ca2+ and IP3 obtained in OpenCOR without the production of IP4. The URL link for the SED-ML le is https://models.physiomeproject.org/workspace/4f9/ le/a779f6ed1fef6592a6646dfd539413b47e56ac04/test_dupont Figure 5 Pathway of interest from the Cooling et al. model [18]. To illustrate the details of the reaction, we show the translocation reaction in this diagram.

Figure 6
Schematic diagram of the F CERIγ signalling pathways The best t for phosphorylation of Grb2. Blue curves show model predictions and red points are observed experimental data in Tsang et al. [23]. The best t for phosphorylation of FCE (left) and phosphorylation of Syk (image not available with this version). Blue curves show model predictions and red points are observed experimental data from Faeder et al. [24].

Figure 9
Fitted model to data from Rajagopalan et al. [25].

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. AdditionalFile1.pdf AdditionalFile2.pdf AdditionalFile3.pdf AdditionalFile4.pdf AdditionalFile5.pdf