The influence of circadian rhythms on CD8 + T cell proliferation upon vaccination: a mathematical modeling perspective

Circadian rhythms have been implicated in the modulation of many physiological processes, including those associated with the immune system. For example, these rhythms influence CD8 + T cell responses that are part of the adaptive immune system. The mechanism underlying this immune-circadian interaction, however, remains unclear, particularly in the context of vaccination. Here, we devise a molecularly-explicit gene regulatory network model of early signaling in na¨ıve CD8 + T cell activation pathway, comprised of three axes (or subsystems) labeled ZAP70, LAT and CD28, to elucidate the molecular details of this immune-circadian mechanism and its relation to vaccination. This is done by coupling the model to a periodic forcing function to identify the molecular players targeted by circadian rhythms, and analyzing how these rhythms subsequently affect CD8 + T cell proliferation dynamics under differing levels of vaccine loads. By performing both bifurcation and parameter sensitivity analyses on the model at the single cell and population levels, we find that applying periodic forcing on molecular targets within the ZAP70 axis is sufficient to creating a day-night discrepancy in CD8 + T cell activation in a manner that is dependent on the bistable switch inherent in CD8 + T cell early signaling. We also demonstrate that the resulting CD8 + T cell activation is dependent on the strength of the periodic coupling as


Introduction
Circadian rhythms, 24 h cycles that modulates many of the body's internal processes (Portaluppi et al., 2012;de Goede et al., 2018), have been linked to the immune system from as far back as the 1960s (Halberg et al., 1960;Esquifino et al., 1996;Kawate et al., 1981;Young et al., 1995).This modulation is thought to be regulated by the presence of molecular clocks, made of clock genes, within cells of the immune system (Houtman et al., 2005;An et al., 2019;De Boer et al., 2001;Coombs et al., 2011;Bollinger et al., 2011;Cuesta et al., 2017;Boivin et al., 2003), including Clock, Bmal1, Per1-3, Cry1-2 ) (Cha et al., 2019).The expression level of these clock genes is regulated by transcription-translation feedback loops that create oscillations in their expression pattern with a ∼24 h period.These clock genes also interact with specific processes within the genome of immune cells and in turn modulate their expression in a cyclical manner, imposing circadian rhythms on them (Haspel et al., 2020;Labrecque and Cermakian, 2015;Jerigova et al., 2022).With evidence suggesting that circadian rhythms affects how the immune system responds to vaccination (Elliott et al., 1972;Long et al., 2016;Phillips et al., 2008;de Bree et al., 2020), it becomes imperative to understand how this relation is manifested molecularly and dynamically at the single cell and population levels.
In innate immunity, the molecular basis of the immune-circadian link is relatively well-established with evidence indicating that these circadian interactions have links to specific processes within innate immune system (Haspel et al., 2020;Labrecque and Cermakian, 2015;Jerigova et al., 2022).Molecular clocks have been characterized in macrophages, neutrophils, eosinophils and natural killer (NK) cells (Keller et al., 2009).Specifically, in macrophages, the clock gene REV-ERBα has been shown to affect the inflammatory function of macrophages by repressing their NF-κB activation pathway (Sato et al., 2014;Wang et al., 2018), whereas in NK cells, mice with Per2 knockdowns displayed reduced levels of proinflammatory cytokines, namely IFN-γ and IL-1β, after an LPS endotoxic shock (Liu et al., 2006).Furthermore, knockouts of Bmal1 in myeloid cells, which make up the majority of innate immune cells, predispose mice to developing acute and chronic inflammation pathologies by disrupting the diurnal variation in monocyte number (Nguyen et al., 2013).
Although significant progress has been made in understanding how innate immunity is influenced by circadian rhythms, the molecular basis of this immune-circadian link remains unclear in adaptive immunity, particularly in CD8 + effector T cell response.While circadian oscillations were observed in the RNA level of lymph node Per2 and circulating naive CD8 + T cells (Fortier et al., 2011;Dimitrov et al., 2009), the molecular source(s) of this diurnic variation remains unknown.A recent transcriptomic study identified an increased activation of certain molecular players within the CD8 + T cell activation pathway such as ZAP70, AKT and mTOR following a daytime vaccination, as well as an increased ability to combat an infectious challenge after a day-time vaccination (both compared to nigh-time vaccination) (Nobis et al., 2019).Moreover, another study has shown an increased level of ZAP70 prior to vaccination in the daytime (Fortier et al., 2011).
These studies thus suggest that certain molecular players within the CD8 + T cell activation signaling pathway interact with circadian rhythms in a way that primes these T cells to exhibit increased activation when vaccinated during the day-time compared to night-time; the details of this interaction, however, remains not fully explored.
In this study, we generate a theoretical framework, through mathematical modeling, to further elucidate the molecular basis of the immune-circadian link in CD8 + T cells following vaccination and its underlying dynamics.We adapted an experimentally validated, CD4 + T cell gene regulatory model (Zheng et al., 2005;Perley et al., 2014) to simulate the immediate signaling events during naïve CD8 + T cell activation following vaccination.
By expanding the model to include a parameter-specific periodic forcing function, we explore the effects of a circadian rhythm on this signaling pathway to determine how it impacts CD8 + T cell proliferation post-vaccination in a molecular-player-specific manner.This is first done via a single-parameter model framework, where bifurcation, sensitivity and time series analyses are performed one parameter at a time to develop a deeper understanding of the dynamics of this modulation.We will then expand this approach to a multiparameter framework to generate more physiologically relevant conclusions about this system.While prior modeling studies have explored the effect of a circadian rhythm on the immune system (Bai and Zhou, 2012;Fan and Wang, 2010;Ji et al., 2010;Wang et al., 2006), these studies were not done using gene regulatory T cell models.
Through our single cell and population level modeling, we demonstrate the crucial role played by the ZAP70 axis (or subsystem) within the immune activation pathway in how a circadian rhythm exerts its effects in modulating CD8 + T cell responses.By first highlighting the presence of a bistable switch in this subsystem, we show that this switch dictates the dynamics of the circadian modulation, and that the latter is dependent on both the vaccine load presented as well as on the magnitude of circadian coupling.Additionally, we establish that this modulation occurs prior to vaccination, suggesting a priming role for a circadian rhythm on CD8 + T cell signaling.Interestingly, our results also reveal that downstream molecular players, such as mTORC1, require their own molecular clocks to produce oscillations, pointing to the existence of a complex, multi-layered circadian influence on CD8 + T cell signaling.Our population modeling further expands on these ideas using a multi-parameter framework to produce more physiological conclusions about the influence of a circadian rhythm on CD8 + T cell activation.

Gene regulatory network model
To explore the dynamics of the circadian rhythm on CD8 + T cell activation, a multi-scale gene regulatory network model of the early events of T cell signaling is considered.The network model is described by a system of ordinary differential equations (ODEs) adapted from literature.This model, first developed in Zheng et al. (2005) and later extended in Perley et al. (2014), describes the immediate signaling events following antigenic stimulation of a single CD4 + T cell.Because of the overlap in these early signaling events between CD8 + and CD4 + cells (Mørch et al., 2020), we have adapted the CD4 4 model to describe CD8 + T cell signaling and modified it to exclude the gene regulatory output signals described in Perley et al. (2014) since CD4 + and CD8 + naïve T cells vary in their gene outputs (Hosking et al., 2014).The model is then coupled to a rhythmic signal representing the circadian rhythm and a single rate-activation parameter is chosen as a molecular target for circadian rhythm coupling.The effects of this coupling on model dynamics are subsequently analyzed.
Briefly, the gene regulatory network model is comprised of the early signaling events following T cell receptor (TCR) phosphorylation (Fig. 1).According to this framework, a phosphorylated TCR (denoted by T CRp) binds to a tyrosine kinase (ZAP 70) and triggers a cascade of interactions that ultimately result in the binding of ZAP 70 to activated Src family protein kinases (SF K).This cascade of interactions is governed by a complex network of feedback mechanisms, regulated by positive feedback loops from the various forms of ZAP 70 (denoted by ZAP * , ZAP p and SF K * -ZAP p), as well as negative feedback loops from Src homology region 2 domain-containing phosphatase-1 (SHP 1), C-terminal Src kinase (CSK), and CREB-binding protein (CBP ).The activated ZAP p and SF K * -ZAP p then feed into a series of reactions governed by various molecular players, such LAT, AKT etc., which ultimately result in the production of mammalian target of rapamycin complex 1 & 2 (mT ORC1 and mT ORC2).Both of these protein complexes, mT ORC1 and mT ORC2, are important in determining T cell fate and effector function (Chi, 2012).The rate equations and parameter values for this model are shown in Appendix 6.1, and a more detailed description can be found in Zheng et al. (2005) and Perley et al. (2014).Furthermore, a framework for this model is shown in Fig. 2 and a list of its parameter values are provided in Table 6.1.

The three key axes of the model
To explore the impact of the circadian rhythm on different components of the model, we divide the gene regulatory network model (Fig. 1) into three subsystems based on how they feed into each other.The first subsystem (labeled the ZAP axis and highlighted by the red box in Fig. 1) includes the reactions that involve ZAP active formation.It begins with the T CRp input, and ends with SF K * −ZAP p and SHP 1.The second subsystem involves the signaling molecules downstream of ZAP active , which include the reactions governing LAT phosphorylation (hence the label LAT axis highlighted by the blue box in Fig. 1).These reactions similarly do not receive feedback from downstream signaling molecules.Finally, the third subsystem involves the signaling molecules included in the CD28 signaling pathway (hence the label CD28 axis highlighted by the red box in Fig. 1)).It starts with the two input signals: ZAP active and CD28, and ends with mT ORC1 * .By exploring the effect of the circadian rhythm on the output of each one of these axes and determining how this circadian influence is filtered from one axis to the next, one can understand which one of them is more sensitive to circadian-induced perturbations, helping us identify potential molecular targets of the circadian rhythm.

Periodic forcing by the circadian rhythm
To investigate the effect(s) of the circadian rhythm on T cell early signaling events, we incorporate a parameter-specific periodic forcing function applied on a given rate constant of one of these early signaling events; that is to say, an external forcing function is imposed onto a parameter of interest within the model, p 0 , to produce a time-dependent variation in its value, denoted by p t .This time-dependent variation in p 0 may get transmitted to downstream signals within the model, producing a time-dependent difference in its output signal.If so, this will allow us to identify which parameter(s) is/are responsible for the day-night discrepancy in CD8 + T cell signaling.
The equation we choose for this periodic forcing is given by where p 0 is the default (original) value of the parameter p that is directly affected by the periodic forcing representing the circadian rhythm, t is time, α p is the parameter-specific 'circadian coefficient' that determines how strongly that parameter is impacted by the oscillation, and β p determines the parameter-specific phase shift of the oscillation.
According to Eqn. 1, the periodic forcing is defined as a sine function that is oscillating with a mean p 0 and a circadian coefficient α that represents the magnitude of circadian rhythm coupling to a specific parameter of interest.
Such description allows us to study the effect of different strengths of circadian coupling on model dynamics.The choice of a sine function to describe the circadian rhythm is consistent with how it is routinely characterized in experimental literature (Elliott et al., 1972;Hickey et al., 1984;Naitoh et al., 1985).
where ZAP active is the percentage of ZAP within a single T cell that has been fully phosphorylated, ZAPp and SFK*-ZAPp are the fully phosphorylated forms of ZAP70, and ZAP total denotes the conserved total quantity of ZAP.
With this output variable, the dynamics of the ZAP axis are explored through parameter sweeps and bifurcation analysis.In Fig. 2(A), we show that for T CRp = 27000 mol, the time course of ZAP active , the output signal of the ZAP axis, changes significantly as the values of the rate activation parameters of the model are altered.For instance, as the rate activation parameter for SHP1, k10f (blue) decreases past a threshold of ∼ 4.7 (mol • min) −1 , the ZAP active output jumps from 0 to a higher steady state of ∼25%.This threshold is also seen with other parameters within the model, as shown with k4f (gray).
To further explore the dynamics of the model in more detail, a oneparameter bifurcation analysis of ZAP active with respect to k10f is conducted, as shown in Fig. 2 Extending this analysis to the LAT and CD28 axes indicates that their dynamics are not as complex.Indeed, parameter sweeps and bifurcation analyses reveal a monostable behaviour across all the parameters within these two subsystems, suggesting that the LAT and CD28 axes do not exert the same level of influence on the overall behaviour of the system as the ZAP axis.This is demonstrated in Fig. 2(C,E), where the time series simulations of P LCγp and mT ORC1 * percentages (relative to their individual total concentration) from the LAT and CD28 axes, respectively, are plotted at different parameter values of k12f and k13f (denoting the forward rates of LAT and P LCγ phosphorylation, respectively) and at different parameter values of k42f and k43f (denoting the forward rates of RhebGDP and mT ORC1 activation, respec-tively) provided in the legends; in these two cases, no threshold-dependent sudden jumps in their steady state output dynamics are observed.This is further confirmed in Fig. 2(D,F) showing a monotonically varying branch of stable equilibria in the bifurcation diagrams of P LCγp and mT ORC1 * precentages with respect to k12f and k42f , respectively.This emphasizes the limited significance of these parameters in defining the dynamics of the LAT and CD28 axes output variables.
This result suggests that a circadian rhythm acting on signaling molecules within the ZAP axis would have the largest impact on model dynamics.This is validated in the next section when performing sensitivity analysis on the parameters of the ZAP axis.To identify sensitive signaling molecules within the ZAP axis, partial rank correlation coefficient (PRCC) of each parameter within a given axis is calculated (Marino et al., 2008).A partial correlation characterizes the relationship between a parameter p i and the output variable, while discounting the effects of other parameters in the model; this is repeated for every activation parameter p within the ZAP axis.This was done using a MATLAB package written by Marino et al. (Marino et al., 2008), which incorporates Latin Hypercube Sampling to generate parameter sets whose partial correlation coefficients are calculated.The results of this are shown in Fig. 3.
From this plot, we designate k4f, k5f, k8f 2, k9f 1 and k10f as sensitive parameters, i.e. with correlation coefficients > |0.6|.This threshold correlation coefficient represents 80% of the correlation coefficient of k10f (whose correlation coeficient is 0.72), where k10f denotes the rate of activation of the SHP 1 molecule.Since SHP 1 has been suggested to possess a circadian component in its regulation (Nobis et al., 2019), we assume that the other circadian-affected parameters of the model should posses a level of sensitivity equivalent or close to that possessed by SHP 1 activation through k10f .
Interestingly, those parameters identified with this approach are involved in the formation of ZAP active , with k4f and k5f denoting the forward and backward rates of SF K * -ZAP p phosphorylation, and k8f 2 and k9f 1 denoting the formation rates of ZAP * and ZAP p, respectively.Additionally, this sensitivity analysis demonstrates the importance of negative feedback in affecting ZAP active through SHP 1 * .Since k10f was the primary negative feedback parameter implicated by this analysis, it will be our target param-eter in the single-parameter analysis conducted in the following section.

The impact of vaccination on the steady state dynamics of the model
Implementing the periodic forcing function defined by Eqn. 1 on k10f (the rate of activation of the SHP 1 molecule) will allow this parameter to exhibit a 24 hr cycle, mimicking the influence of the circadian rhythm.This is shown in Fig. 4(A), where k10f t is plotted according to the circadian time of day, with [0, 12] hr representing day-time and [12,24] hr representing night-time.
The amplitude of this diurnal variation in k10f t varies according to α k10f denoting the coefficient of variation; when α k10f = 0.5 (solid line), larger oscillations are produced compared to α k10f = 0.25 (dashed line).The respective day and night time maxima and minima of k10f t are identified by an orange (yellow) and blue (light blue) dot for each value of α k10f , respectively.Recall that (see Fig. 2) the bifurcation diagram of the model exhibits a bistable switch with two saddle-nodes (SN1 and SN2) and two Hopf bifurcations (HB1 and HB2) at which two envelopes of unstable periodic orbits emerge and eventually terminate at homoclinic bifurcations (HM1 and HM2, respectively).Superimposing the orange and blue dots representing the maxima and minima of k10f t when α k10f = 0.5 on this bistable switch, as shown in Fig. 4(B), the orange dot ends up lying on the upper branch of stable equilibria (to the left of HB1), while the blue one ends up lying on the lower branch (to the right of HB2).This means that during a 24 hr day-night cycle in k10f t , the solution trajectory of ZAP active will oscillate between these two dots that represent the extrema of k10f t and cause it to jump back and forth between the two stable branches of equilibria.The jumping to the lower branch will occur during the transition from daytime to nighttime when the trajectory moves rightward along the upper branch of ZAP active stable equilibria until it crosses HB1, leading the trajectory to jump to the lower stable branch and causing ZAP active to decrease.Conversely, during the transition from nighttime to daytime, the trajectory will move leftward along the lower branch of ZAP active stable equilibria until it crosses HB2, leading the trajectory to jump to the upper stable branch and causing ZAP active to increase.
Thus, the circadian rhythm on ZAP axis activation parameters can have a drastic effect on ZAP active output.This drastic change in dynamics occurs because of the presence of the bistable switch with its two Hopf bifurcations: HB1 & HB2, lying within the parameter space afforded by α k10f = 0.5.This cycling between the two branches of equilibria every 12 hr suggests that the circadian rhythm can cause periodic oscillations in ZAP active with large am-plitudes.This then posses the question of how such mechanism interacts with vaccines?
To answer this question, we will first investigate how the impact of the circadian rhythm on T cell activation is affected by the T CRp (the number of phosphorylated TCRs), which we assume to be an indicator of vaccine strength.This is done by considering two other values of T CRp: one that is higher than the one used in Fig. 4 According to this realization, solution trajectories are unable to cross the Hopf bifurcation HB2 and solely oscillate along this lower branch with very low amplitude, i.e., without exhibiting drastic day-night shift in ZAP active .
Such realization thus indicates that if the influence of the circadian rhythm is not strong enough, then an intermediate level of vaccination may still not be able to generate large amplitude day-night cycle in ZAP active .
Taken together, these results thus demonstrate that the day-night discrepancy in signaling is dependent on both having a sufficiently large T CRp input, as well as on possessing a sufficiently high α k10f , such that the k10f t is able to cross HB1 and HB2 of the model.This further highlights the complex relationship between vaccine load and the strength of circadian modulation when mediating CD8 + T cell response upon vaccination.Further experimental work is needed to validate the extent of this relationship.To address this latter issue, we hypothesize that upon vaccine application the effect of the circadian oscillations on k10f is significantly attenuated with a much lower value for α k10f .To test this hypothesis, we set α k10f = 0 right when vaccination is applied at both daytime (at t = 12 hr) and nighttime (at t = 24 hr) as shown in Fig. 5(D).This is done by freezing (i.e., fixing) the value of k10f at its respective value upon vaccine onset.Even though vaccine application still transitions the system from the state S 0 to the state S 1 , repetitive large amplitude oscillations no longer occur -only small amplitude oscillations in ZAP active is seen prior to vaccination due to the circadian clock as described before for S 0 .Interestingly, nighttime vaccination causes a very moderate increase in ZAP active (Fig. 5(E)), whereas daytime vaccination causes a significant increase in ZAP active (Fig. 5(F)), reaching a significantly elevated steady state due to the transition to the upper branch of stable equilibria highlighted in Fig. 5(B).The lack of oscillations in k10f t post vaccination prevents the solution trajectory from crossing HB1 and returning to the lower branch.These results thus corroborates the experimental finding that daytime vaccination is more effective in triggering T cell activation, and suggests that for this to be true, the diurnal effect of the circadian rhythm on early events of T cell signaling is significantly attenuated post vaccination.
Taken together, these results indicate that the circadian rhythm primes or facilitates the activation of naïve CD8 + T cells to produce the day-night discrepancy in ZAP active output, rather than directly controlling it.

Characterizing the filtering properties of the ZAP axis
Throughout this study, we have explored the implications of the circadian rhythm on the ZAP axis, which we have found to be the primary driver of model dynamics with its intrinsic bistable switch.Here, we expand this study by exploring how the effects of the circadian rhythm acting directly on the ZAP axis is transmitted to the downstream LAT and CD28 axes within the model in the absence of vaccination.Doing so would indicate whether a circadian rhythm within the ZAP axis is capable of producing a day-night discrepancy in the signaling of other downstream molecular players associated with these two axes.If not, this may suggest the presence of additional circadian clocks acting directly on these downstream signaling molecules.
To explore how the influence of the circadian oscillation within the ZAP axis in the absence of vaccination is transmitted to other axes of the model, we now consider the expanded model of Fig. 2 that includes the LAT and CD28 axes.In order to quantify how the circadian-dependent oscillations of k10f t are transmitted to other variables in these two axes , we will use the 'oscillatory contribution' ratio, ω i , for every variable (molecular player) i of the model.This ω i is defined as the relative effect of the oscillation in k10f t to the mean of the variable i at steady state.It is given by where i represents the variable of interest, SS amp,i represents the amplitude of the steady state oscillations, and SS mean,i represents the mean of the variable.A larger ω i indicates that a larger proportion of a variable's steady state can be explained by k10f t , whereas a smaller ω i indicates that the circadian-dependent oscillations in k10f t minimally alter the steady state of the variable i.To quantify ω i , we will always assume that k10f t is oscillatory with a circadian coefficient of α k10f = 0.5.
By plotting the magnitude of ω i for the major variables within the ZAP axis, as shown in Fig. 6(blue bar), we find that SHP 1 active is associated with the largest ω SHP 1 , while the other variables (including i = SF K * , CSK * , CBP p) exhibit a reduced ω i .This is expected since a periodic forcing on k10f directly regulates SHP1 activation, altering SHP 1 * .
For the variables within the LAT axis, on the other hand, Fig. 6(B) shows that the magnitude of ω i (i = LAT p, P LCγp, SOSb) are lower than those for the key molecular players within the ZAP axis: ZAP active and SHP 1 * , suggesting that some of the k10f -induced oscillation is filtered out between axes.
However, the oscillations remain identifiable within the axis, indicating that a circadian rhythm in the ZAP axis is sufficiently transmitted to the LAT axis.
Interestingly, in the CD28 axis, Fig. 6(C) shows that while the oscillation from k10f t is adequately transmitted to many of the signaling molecules (including i = P I3K * , P DK1 * , AKT * , mT ORC2 * ) with moderate values of ω i , it is completely filtered out when i = T SC1 − 2, RhebGT B, mT ORC1 * (i.e., ω i = 0 for all of these cases).This suggests that oscillations in k10f t directly influencing the ZAP axis may not be able to produce the experimentally observed oscillations within mT ORC1 * and its immediate precursors (Nobis et al., 2019), despite producing them in the other signaling molecules of this axis.
This was resolved by imposing a periodic forcing component, representing the influence of the circadian rhythm, directly on RhebGT P activation parameter k42f with α k42f = 0.5; this resulted in a high ω RhebGT P and ω mT ORC1 active of 49% and 49.5%, respectively, as shown in Fig. 6(C, inset).
Curiously, adding a periodic forcing onto any of the more upstream signaling molecules within the CD28 axis, such as AKT, PI3K or even TSC1-2, was not able to produce oscillations in RhebGT P and mT ORC1 * , further indicating the presence of an independent circadian regulation of RhebGT P and mT ORC1 * .Since S6 (a molecular derivative of mT ORC1), has been experimentally shown to have a circadian-influenced day-night discrepancy (Nobis et al., 2019), these modeling results suggest that a circadian rhythm within the ZAP axis may not be sufficient to explain these downstream oscillations, pointing to a potential secondary molecular clock acting directly on the gene transcripts of mT ORC1-related players.

Population model framework
To better understand T cell responses to circadian oscillations and vaccines, it is imperative to obtain a broader view of model dynamics at the population level.We do so by simulating the steady state dynamics of k = 100 model realizations of the ZAP axis, each with a randomly selected parameterization of the model.This is done after eliminating the oscillatory solutions using a method called parallel coordinate plot (see Appendix 6.3).
We initially conduct this parameterization by first targeting one parameter, namely, k10f k .This is done by randomly assigning random values for this parameter within the range [0.5×k10f 0 , 1.We refer to this plot as a 'macro-bifurcation' diagram.
According to this macro-bifurcation diagram, we can split the parameter space of the model into 3 regions labeled insensitive, sensitive and elevated.
This distinction depends on the distribution of ZAP active at each value T CRp.
The insensitive region is the range of T CRp ∈ [1, 2]•10 4 mol, where ZAP active is purely in its inactive state, regardless of the parameterization of k10f .The sensitive region, on the other hand, corresponds to the range of T CRp ∈ [2, ∼ 3.10] • 10 4 mol, where ZAP active can be either elevated or not, depending on the value of the parameter k10f t .We denote this region as possessing 'macro-bistability', as the magnitude of its output depends on the value of the circadian-affected parameter k10f .Finally, the elevated region is the one that corresponds to the range T CRp ∈ [∼ 3.10, 3.5] • 10 4 mol, where ZAP active is purely elevated, regardless of the value of k10f t .
Unique to the macro-bifurcation diagram is the ability to expand this analysis to n parameters, whereas our previous analysis was limited to 1 or 2 dimensions using the one and two-parameter bifurcations.Indeed, by expanding this approach to include variability in the positive feedback parameters k4f, k9f 1, negative feedback parameters k7f, k10f, k11f and our manually derived threshold parameter k8f 1 ((see Appendix 6.3.2)) as shown in Fig. 7(B), we can obtain macro-bifurcation diagram of ZAP active with respect to T CRp ∈ [0.5, 5.5] • 10 4 mol.We have chosen these parameters because they can act as possible molecular targets of the circadian rhythm in early CD8 + T cell signaling and because they were a subset of 'sensitive' parameters identified by our sensitivity analysis of Fig. 3 (for the full definition of k4f, k9f 1, k7f, k10f, k11f, k8f 1, please see Tables 6.1).Combining their effects in a macro-bifurcation diagram can thus provide a more physiological understanding of how the circadian rhythm can impact CD8 + T cell activation.Simulating 100 realizations of the model, by randomly sampling these parameters from their respective ranges defined by [0.5 × par 0 , 1.5 × par 0 ], where par = k4f, k9f 1, k7f, k10f, k11f, k8f 1 and par 0 is the default parameter value of each one of them (see Tables 6.1 T cell signaling in response to vaccines.Further analysis, coupled with experimental work, may help determine how additional rate activation parameters can alter the macro-bifurcation diagram, providing a more targeted approach to explore the effects of the circadian rhythm on CD8 + T cell activation.

Discussion
In this study, we provided detailed analysis of the dynamics of the signaling molecules involved in producing the day-night discrepancy seen in early CD8 + T cell signaling.Additionally, we explored how this day-night discrepancy may be altered by vaccines of different strengths.This was done through a mathematical modeling approach, which provided flexibility in studying how changes in rate activation parameters can affect CD8 + T cell activation.
More specifically, a reference model for CD4 + T cell activation obtained from Zheng et al. (2005); Perley et al. (2014) was adapted to CD8 + T cell activation and used to conduct our analysis.This model was then split into 3 subsystems, and a dynamical systems analysis was conducted on each one to identify the main drivers of model dynamics.Through this analysis, the subsystem associated with ZAP 70 activation (ZAP axis) was shown to drastically affect the dynamics of the other two downstream LAT and CD28 axes by uniquely possessing a bistable switch.This is in agreement with prior modeling literature suggesting a bistability within ZAP70 signaling (Lipniacki et al., 2008).To further pinpoint potential molecular targets of the circadian rhythm within the ZAP axis, a global sensitivity analysis was conducted using Latin hypercube sampling and partial rank correlation coefficient (Marino et al., 2008).This was done to identify how parameter perturbations within the ZAP axis can affected the dynamics of its output signal, namely ZAP active .Since the circadian rhythm was implemented in a parameter-specific manner, identifying the most sensitive parameters of the model provided insights into potential targets of the circadian clock in CD8 + T cell signaling.Through this analysis, we identified that ZAP 70 and SHP 1 activation parameters, namely k4f, k5f, k10f, k9f 1, k8f 2 were the most sensitive to perturbations, conclusions that are aligned with the results of transcriptomic analyses of CD8 + T cells in mouse genomes (Nobis et al., 2019).
To analyze the effects of this parameter-specific circadian coupling on model dynamics, a single-parameter bifurcation analysis approach was applied.Due to the bistable switch within the ZAP axis, we showed that, for a given vaccine input T CRp, an elevated ZAP active was dependent on the value of the circadian-coupled bifurcation parameter k10f t , especially relative to its Hopf bifurcation points, thus producing a day-night discrepancy.
This day-night discrepancy was also shown to be affected by the vaccine load presented to the system, via the input parameter T CRp.This suggests a dependency between the vaccine load presented to CD8 + T cells and the influence of the circadian rhythm, in such a way that there exists a threshold in T CRp below of which the circadian rhythm may not affect CD8 + T cell signaling.These results suggest a complex relationship between vaccine load and circadian coupling, which may explain why characterizing the vaccination effects of circadian rhythm in T cells has been relatively difficult thus far.
Through a time-series analysis on the single-cell scale, we showed that only a circadian rhythm in the S 0 state of the model produces a physiologically valid day-night discrepancy and that such circadian effect must abort post vaccination.This suggests that the physiological circadian rhythm may act during the pre-vaccination state of a CD8 + T cell, 'priming' them to respond to a vaccine.This result corroborates previously observed experimental results showing that CD8 + T cells are more primed for activation even prior to vaccination (Nobis et al., 2019).Additionally, this finding is aligned with research in macrophages showing that the Bmal1 clock gene expression is reduced after an immune challenge via LPS (Curtis et al., 2015).
To expand this analysis to other molecular players of CD8 + T cell signaling, the filtering properties of molecular players within the ZAP axis to the downstream LAT and CD28 axes were explored.Our results re-vealed that an S 0 circadian rhythm within the ZAP axis is not transmitted to certain downstream signaling molecules within the CD28 axis, namely T SC1 − 2, RhebGT P & mT ORC1.However, a circadian rhythmicity to their activity was restored by adding a circadian component to their own activation parameters, namely to k42f .Since a day-night discrepancy was experimentally observed in mT ORC1 * (Nobis et al., 2019), our results suggest that this day-night discrepancy may be due to a secondary, independent molecular clock acting directly on mT ORC1, rather than being due to the activity of the upstream molecular clock within the ZAP axis.
To produce a more physiological representation of CD8 + T cell responses to vaccination and how that is impacted by circadian rhythm, we implemented a multi-dimensional population level modeling approach.This allowed us to explore the effect of the circadian rhythm on multiple molecular players simultaneously, as well as determine how heterogeneity in CD8 + T cell population affect their robustness to vaccination.Through these simulations, three regions of activity with respect to vaccine load defined by T CRp were identified when key circadian-dependent parameters were randomly sampled within their respective ranges: insensitive, sensitive (macro-bistable) and elevated.This macro-bifurcation diagram may provide a framework to explore the results of adding a circadian perturbation to additional molecular players with the purpose of motivating further experimental research.
One limitation of this study is its assumption that changes in activated  (2001,2003); Coombs et al. (2002); François et al. (2013).That said, this study provided several novel insights into the effect of the circadian rhythm on CD8 + T cell signaling, in the context of vaccination, and serves a valid framework for assessing the influence of the circadian rhythm in other models.

Model Parameterization
Here, we discuss the key reactions and framework of the model as well as note any modifications made by us to accommodate CD8 + T cell signaling.
Broadly. the model encompasses the immediate events following antigen presentation, on a single T cell scale.The equations of the model are shown in Table 2.

Variable Differential Equation
ZAP axis CD28 axis Table 1: Differential equations of the various signaling molecules within the model.Equations adapted from Zheng et al. (2005) and Perley et al. (2014).

Components
ZAP axis     invalid result of the model, as this behaviour is unlikely to exist.To resolve this issue, we have developed a systematic framework to eliminate this oscillatory regime while maintaining the other model regimes.This is the focus of the following section.

Eliminating the intrinsic oscillatory regime
This approach is adapted from the Dynamics Visualisation based on Parallel Coordinates (DYVIPAC) framework by Nguyen et al. (2015).Broadly, it involves a population-based approach to uncover system's dynamics.The m parameters whose effects on model dynamics we are interested in are first selected, and a range of values for each parameter is set.A sample of n model realizations (i.e.simulated steady state response) is then created, each with a randomly selected parameterization from within the pre-determined ranges.These realizations are integrated starting from the same initial conditions, and their steady state solution is classified according to its dynamics.
A parallel coordinates plot is then created linking each realization of the model to its respective steady state classification.This allows us to obtain a multi-dimensional view of model dynamics, and allows us to reliably identify parameter regimes where the oscillatory solution exists.By limiting our parameter sampling ranges to exclude such parameter regimes, the oscillatory solutions can be excluded from our model.This approach is explored in more detail in the following sections.

Network model and parameter selection
The first step of this framework is to identify the parameters of interest, or more specifically those that we vary in this study.These parameters are T CRp and k10f ; k8f 1 is also selected since it allows us to more easily eliminate the oscillatory regime.The number of parameters defines the dimensionality m of this approach (i.e., m = 3).

Parameter sampling
With the dimensionality m of the framework determined through parameter selection, an appropriate sampling range of these parameters now needs to be set.Each realization of the model will have a randomly selected value, selected from within this range.Since we are interested in the effect of the circadian rhythm, this parameter range should include the degree of modulation due to the circadian rhythm.This range is therefore defined by where p i is the the value of a parameter i randomly sampled from within a set range.This p i , when sampled over a sufficiently large k realizations, demonstrates all the possible values this parameter can have when subjected to a circadian rhythm.p i,0 denotes the default parameter value, and p i,min and p i,max denote the minimum and maximum endpoints from which this parameter can be sampled from.

Setting a physiological p i range
To impose a physiologically-relevant limit for the influence of the circadian rhythm, we note that in CD8 + T cell signaling it was shown that IFNγ, along with other CD8 + T cell cytokines, approximately vary within a [0.5×p 0 , 1.5× p 0 ] range depending on the time of initial vaccination (Fortier et al., 2011;Nobis et al., 2019).We will therefore use this range to describe the influence of the circadian rhythm on CD8 + T cell signaling.To produce this p range, we choose the range of viable circadian coefficient α to be α ∈ [0, 0.5].It should be noted that this only defines the maximum range of p; some realizations of the model may have a parameter subjected to a weaker circadian influence (defined by an α < 0.5).
It is also important to note that the above range assumes that a [0.5×p 0 , 1.5× p 0 ] discrepancy in cytokine production can be traced linearly to a discrepancy of a similar magnitude in ZAP signaling.While we show in this study that the axes downstream of ZAP do not contribute significantly to model dynamics, we acknowledge that there may be additional factors contributing to the discrepancy experimentally seen in IFNγ and other cytokines.Regardless, this assumption is made in this study to facilitate our analysis and selection of a viable parameter range.

Analysis of dynamical behaviour
With k model realizations selected from parameters within physiological ranges, the individual realizations can be integrated and their steady-state solution analyzed and classified according to its properties.Whereas DYVIPAC uses model eigenvalues to calculate stability (Nguyen et al., 2015), our approach classifies solutions based on their steady state behaviour.Namely, we classify solutions as either elevated (i.e., corresponding to those with a higher than 5% active ZAP concentration), oscillatory (i.e., corresponding to those with a slope greater than a manually selected threshold of 0.0053) or inactive (i.e., corresponding to those fulfilling neither of the previous criteria).
We have found that performing the classification in this manner is effective due to two important dynamic features of the model, including: (i) its fast dynamics, allowing it to reach steady state relatively quickly, and (ii) the relaxation oscillation-like behaviour it exhibits within the oscillatory regime with its rapid up and down-stroke oscillations that facilitate our ability to classify it according to its slope.
With this, we are left with a matrix of k model realizations, each with its unique paramaterization of the m circadian-affected parameters, and with a label designating the classification of each k i solution based on the properties of its steady-state dynamics.

Parallel coordinate plot visualization
This dataset is visualized using a parallel coordinate plot, which displays how these different parameter values can influence the dynamics of the model.The parameter values are scaled between [0, 1] according to their respective maximum and minimum sample values.This normalization is generated using the equation where, p i,norm is the normalized parameter value between [0, 1] of a specific parameter of interest.p i is the parameter value of a particular model realization, and p i,max and p i,min are the respective maximum and minimum parameter values across all n realizations.This visualization allows us to obtain a broader picture of model dynamics and how these dynamics can be shaped by different parameter combinations.

Eliminating the oscillatory regime
Using the parallel coordinate plot, the oscillatory regime can be extracted from the original scaled range of an indicator parameter p i by removing the continuous interval within this range associated with the oscillatory behaviour (i.e., the interval that has segmented yellow lines).The indicator parameter p i is then re-scaled according to the remaining interval and sampling is then restricted to this constrained range (which by design excludes oscillatory solutions).This is done in Fig. S2(B), where the elevated and inactive solutions from the previously presented parallel coordinate plot are hidden, to better visualize the oscillatory regimes.
We identify that oscillatory solutions seem to be confined to a normalized k8f 1 range of [0.23, 1] of k8f 1.Thus, by only sampling from k8f 1 ∈ [0, 0.23], the oscillatory regime can be completely avoided.
Using this threshold of k8f 1 max = 0.23 (mol • s) −1 , we reverse the previously described normalization process, to obtain the raw maximum value of k8f 1 max = 6.12 (mol • s) −1 .To obtain the new p i,0 for k8f 1, the trivial calculation of p i,0 = k8f 1 max /1.5 is performed.The preceding calculation is repeated to find k8f 1 min , thus producing a new range k8f 1 ∈

Figure 1 :
Figure 1: Schematic of the gene regulatory network adapted from Zheng et al. (2005) and Perley et al. (2014), and considered in our model framework to describe early signaling in CD8 + T cells.The model can be divided into three main subsystems: the ZAP axis (red box), LAT axis (blue box) and CD28 axis (green box).Notice how the LAT and CD28 axes are linked to the ZAP axis in a feedforward manner without any feedback, and that the LAT and CD28 axes are similarly not linked in a feedforward manner.Asterisk is used to represent the active forms of the relevant molecules.
3.1.2.Parameters of ZAP 70 and SHP 1 activation are the most sensitive to perturbation Fig. 4(B).Indeed, when T CRp = 40000 mol, the two dots become part of the upper branch of stable equilibria, producing an elevated output of ZAP active independent of time of day, whereas when T CRp = 10000 mol, the inverse is seen, with the two dots lying on the lower branch of stable equilibria, producing baseline output of ZAP active independent of time of day.Thus there seems to be an intermediate range of T CRp within which the circadian rhythm can produce the largest day-night discrepancy in ZAP active .This relationship between T CRp and k10f is further demonstrated by plotting the two-parameter bifurcation of ZAP active with respect to both T CRp and k10f as shown in Fig.4(E).The region bounded by the Hopf bifurcations (blue lines)demarcate the bistable regime, whereas the upper and lower regions denote the monostable elevated and baseline states of the model, respectively.Notice how the bistable regime shifts horizontally as the value ofT CRp is varied, indicating that there is an intermediate range of T CRp within which the largest day-night discrepancy in ZAP active can be generated.This means that the level of vaccines administered can alter T CRp and cause the regime of bistability to shift accordingly.Now that we have explored how vaccines may affect bistability with respect to T CRp and k10f , we now examine how changing the the strength of the circadian rhythm through α k10f can affect outcomes.By taking T RCp = 27000 mol and setting α k10f = 0.25 to produce the two extrema of k10f t identified by the yellow and light blue dots in Fig.4(A), we find that superimposing these dots on the bifurcation diagram of ZAP active with respect to α k10f (see Fig.4(F)) places them on the lower stable branch of equilibria.
3.3.Impact of circadian coupling on temporal evolution of CD8 + T cell signalingNow that we have determined the underlying steady state dynamics of circadian modulation of T cell activation upon vaccination, we now explore how the time series of the ZAP axis is affected by k10f t and the potential physiological implications of this on CD8 + T cell signaling.To do so in the context of vaccination, we define two possible states that can be attained by the model: S 0 and S 1 .S 0 represents the starting state of the model (prevaccination) and corresponds to a given initial amount of T CRp, assumed to be T CRp 0 = 100 mol; a nonzero basal level of T RCp is chosen for this state to accommodate pMHC-TCR contacts with self or non-cognate antigens.S 1 , on the other hand, denotes the post-vaccination state and corresponds to when T CRp >> 10000 mol (>> T CRp 0 ).The effect of the circadian rhythm on the time series of the ZAP axis in each of these states is then analyzed (see Fig.5) by exploring the effects of a k10f -induced circadian rhythm on both S 0 or S 1 .Imposing a 24-hr circadian cycle on k10f in the form of a sine wave with α k10f = 0.5 as in Fig.5(A) and initiating model simulations from S 0 induce a repetitive very small-amplitude oscillations is ZAP active (Fig.5(B, inset)).Although T RCp = 10000 mol in Fig.5(D), the monostable dynamics detected in this figure is responsible for generating these small amplitude oscillations that shuttle ZAP active between the orange and blue dots.This is because the bifurcation diagram at S 0 is essentially the same as the one in Fig.5(D).However, when vaccination is applied, the state of the system switches to S 1 .By setting T CRp = 29000 mol in S 1 , ZAP active then produces repetitive large amplitude oscillations regardless of the time of vaccination (Fig.5(B,C)), with the daytime vaccination (at t = 12 hr) producing a faster response compared to the nighttime vaccination (at t = 24 hr).These large amplitude oscillations are reminiscent to those cycles between the orange and blue dots highlighted previously in Fig.5(B); solution trajectories, in this case, cross the two Hopf bifurcations back and forth and jump between the branches of stable equilibria.Given that such oscillations indicate repetitive activation and deactivation of T cells, then this would preclude them from being physiological.
5×k10f 0 ] for each model realization k, where k10f 0 denotes the default parameter value of k10f listed in Table6.1.Doing so generates heterogeneity between each model realization through this parameter.These model realizations are then all integrated from the same initial conditions at increasing values of T RCp ∈ [1, 3.5] • 10 4 mol and the steady state values of ZAP active are plotted with respect to T RCp (along the x-axis) as shown in Fig.7(A).Each point on the graph denotes the steady state ZAP active at a certain parameter value of k10f k within the previously defined range.This provides a snapshot of 100 steady states of the model, each at a different circadian time defined by the value of k10f k ∈ [0.5 × k10f 0 , 1.5 × k10f 0 ] (recall that we assumed the circadian rhythm affects the value of k10f when we defined k10f t ).Plotting the steady states of these 100 realizations with respect to T CRp, the surrogate for vaccine strength, highlights how the model output varies under differing vaccine loads.This produces what looks like a one-parameter bifurcation diagram that encompasses not only the effects of the circadian rhythm on T cell activation through k10f , but also the effects of vaccine strength through T RCp.
Fig. 7(A): T CRp ∈ [2, 3] • 10 4 mol in panel (A) vs T CRp ∈ [2, 4] • 10 4 in panel(B).This demonstrates that our theoretical framework is maintained when adding additional parameters, and can provide insights into how adding additional parameters to our macro-bifurcation analysis may affect early CD8 + ZAP70 and mTOR translate to changes in CD8 + T cell proliferation.Our future work involves expanding this model to cover the downstream CD8 + gene regulatory outputs.This would expand the scope of the model and allow us to make direct predictions about cytokine production and subsequent proliferation of CD8 + T cells.This expansion would also assist in the experimental validation of the model, since outcomes would be closely linked to CD8 + T cell responses.Alternatively, a different reference model may be chosen as a starting point for this analysis using previously proposed mathematical models of T cell signalingLipniacki et al. (2008); De Boer et al.

6. 2 .
Defining a valid parameter space for the modelWith the parameterization of the model described inZheng et al. (2005) andPerley et al. (2014), the model is able to produce intrinsic oscillatory signals, shown in Fig.S1.For instance, setting T CRp = 20000 mol and k7f = 4.89 • 10 −7 (mol • s) −1 and plotting the one-parameter bifurcation of ZAP p with respect to k10f (the rate constant for SHP 1 activation), as shown in Fig.S1(A), produces a diagram that effectively eliminates the bistable behavior of the model and produces an oscillatory regime defined by the envelopes of stable limit cycles initiated from the saddle-node of periodic orbits (SNP) and terminate at a saddle-node on an invariant circle (SNIC) bifurcation (formed by SN2 and a homoclinic bifurcation HM).The stable envelopes of limit cycles produce periodic solutions whose time series simulations exhibit periodic spiking in ZAP active , as shown in Fig.S1(B).This spiking pattern in ZAP active has not been observed experimentally, indicating that this is a non-viable parameter regime of the model.Phys-iologically, since ZAP 70 remains bound to the TCR throughout antigen presentation, this spiking pattern suggests that ZAP 70 is repeatedly getting phosphorylated and subsequently dephosphorylated in the time scale of seconds, all while remaining bound to the phosphorylated T CR (via the non-zero T CRp parameter).This result contradicts the literature showing a continuous level of phosphorylated ZAP 70 throughout antigen stimulation Houtman et al. (2005); An et al. (2019).The lack of evidence showing a spiking ZAP 70 pattern as well as the contrast to experimentally observed ZAP 70 binding/unbinding kinetics leads us to designate this regime as an Each model realization is linked to its respective steady state dynamics according to the value of its parameters of interest.This visualization is shown in FigureS2(A), where a multi-dimensional analysis of three parameters, including T CRp, k10f (which represents SHP1 activation) and k8f 1 (which represents ZAP activation), is conducted.In this analysis, k = 1000 model realizations are used, where each realization is represented individually by a colored segmented line linking the parameter values generating this realization.The parameter values of each realization lie above the ticks on the x-axis representing each parameter, and are connected by a segmented line color-coded according to the stability property of that particular model realization (i.e., the stability of the steady-state solution associated with that realization).Three colors are used, blue, red and yellow representing inactive, elevated and oscillatory solutions, respectively.

Figure 3 :
Figure3: Partial rank correlation coefficients (PRCC) of the activation parameters of the ZAP axis when correlated to ZAP active output.A higher PRCC value indicates a higher correlation between values of that parameter with ZAP active (positively correlated if PRCC > 0, negatively correlated if PRCC < 0).The following parameters k4f, k5f, k8f 2, k9f 1 and k10f are identified as sensitive with correlation coefficients > |0.6|.

Figure 4 :
Figure 4: (A) Time series plots of k10f t when α k10f = 0.5 (solid line) and α k10f = 0.25 (dashed line) showing 24 hr cycles with different amplitudes.Orange (yellow) and blue (light blue) dots show the diurnal maxima and minima of k10f t when α k10f = 0.5 (α k10f = 0.25).(B) One-parameter bifurcation of ZAP active with respect to k10f at T CRp = 27000 mol.Details of the bifurcation diagram are described in Fig. 2. Orange and blue superimposed dots represent the respective day and night extrema of k10f t identified in (A).Note the abrupt change in steady state upon crossing HB1 and HB2.(C), (D) Oneparameter bifurcation of ZAP active with T CRp = 40000 mol (C) and T CRp = 10000 mol (D).Note the monostable response across the k10f t parameter space.(E) Two-parameter bifurcation of ZAP active with respect to k10f and T CRp.Red lines denote continuations of saddle-node bifurcations (SN1 and SN2), while blue lines denote continuations of the Hopf bifurcations (HB1 and HB2).Bistability is bounded by HB1 and HB2.(F) one-parameter bifurcation of ZAP active at T CRp = 27000 mol as in (B), but with α k10f = 0.25; the day-night maxima identified in (A) are superimposed as yellow and light blue dots.

Figure 5 :
Figure 5: (A) Time series plot of k10f t when α k10f = 0.5 showing 24 hr cycles.Dasheddotted lines indicate the moments of vaccine applications at daytime (t = 12 hr) and at nighttime (t = 24 hr).(B) Time series simulations of ZAP active prior and post vaccination applied at nighttime (t = 24 hr).Inset: magnification of the time series simulations prior to vaccination to highlight the small amplitude oscillations associated with the state S 0 .(C)Time series simulations of ZAP active prior and post vaccination applied at daytime (t = 12 hr).Notice the repetitive large amplitude oscillations associated with the state S 1 produced post vaccination at nighttime and daytime; the onset of the large amplitude oscillations in the former case (nighttime vaccination) is delayed relative to the latter case (daytime vaccination).(D) Time series plot of k10f t when α k10f = 0.5 showing 24 hr cycles that get interrupted upon vaccination at daytime (orange; 12 hr) and nighttime (blue; 24 hr).When vaccination is applied in each case, k10f t is then kept steady at its respective value.(E) Time series simulations of ZAP active prior and post vaccination applied at nighttime (t = 24 hr).(C) Time series simulations of ZAP active prior and post vaccination applied at daytime (t = 12 hr).Inset: magnification of the time series simulations prior to vaccination to highlight the small amplitude oscillations associated with the state S 0 .Notice the much larger response associated with S 1 produced post vaccination at daytime compared to nighttime, as well as the absence of oscillations post vaccination in both cases.

Figure 6 :
Figure6: Bar plots of the magnitude of the oscillatory contribution ω i for different molecular players within the (A) ZAP axis, (B) LAT axis, and (C) CD28 axis when k10f t is oscillatory with α k10f = 0.5.The larger the ω i value, the stronger the influence of the circadian-induced oscillations have on the steady state level of the molecular player i of interest.Inset in panel (C): Bar plots of the magnitude of the oscillatory contribution ω i for the molecular players i = RhebGT P and mT ORC1 * within the CD28 axis when k42f t is oscillatory with α k42f = 0.5.

Figure S2 :
Figure S2: Parallel coordinate plot of T CRp, k8f 1, k10f .(A) The three parameters are displayed on the x-axis and their normalized parameter values are displayed on the y-axis.Each segmented colored line represents a given model realization (steady state response) sampled randomly from within a pre-determined range.The segmented lines are colorcoded according to the type of steady state obtained: blue for inactive, red for elevated and yellow for oscillatory solution.(B) The isolated oscillatory solutions (yellow) from (A). (C) The analysis repeated with a constrained parameter range for k8f 1.Notice the lack of yellow oscillatory solutions.

Table 2 :
Perley et al. (2014)f the various signaling molecules within the model.Equations in ZAP axis are multiplied by 60 to account for time-scale differences in parameter values.Model equations adapted fromZheng et al. (2005)andPerley et al. (2014).