The influence of circadian rhythms on CD8+ T cell activation 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 within 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 the näı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 activation under differing levels of T cell receptor (TCR) phosphorylation, which we designate as vaccine load. 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 create 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 well as on the level of TCR phosphorylation. Our results show that this day-night discrepancy is not transmitted to certain downstream molecules within the LAT subsystem, such as mTORC1, suggesting a secondary, independent circadian regulation on that protein complex. We also corroborate experimental results by showing that the circadian regulation of CD8+ T cell primarily acts at a baseline, pre-vaccination state, playing a facilitating role in priming CD8+ T cells to vaccine inputs according to time of day. By applying a population level analysis using bifurcation theory and by including several hypothesized molecular targets of this circadian rhythm, we further demonstrate an increased variability between CD8+ T cells (due to heterogeneity) induced by its circadian regulation, which may allow a population of CD8+ T cells to activate at a lower vaccine load, improving its sensitivity. This modeling study thus provides insights into the immune targets of the circadian clock, and proposes an interaction between vaccine load and the influence of circadian rhythms on CD8+ T cell activation. Highlights Potential targets of circadian rhythms within the ZAP70 signaling pathway were identified. The level of vaccine load to a CD8+ T cell was shown to be crucial in dictating a circadian rhythm’s influence on its signaling response. The ’priming’ effect of a circadian rhythm on CD8+ T cell activation upon vaccination was demonstrated. mTORC1 and its immediate upstream signaling molecules were shown to be regulated by circadian rhythms through independent mechanisms. A heterogeneity in CD8+ T cells, induced by its circadian regulation, may influence their sensitivity to vaccination.


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 ) (Dibner et al., 2010).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 modulate the expression of specific molecular targets within immune cells in a cyclical manner, imposing circadian rhythms on them (Labrecque and Cermakian, 2015;Jerigova et al., 2022).With evidence suggesting that circadian rhythms affect how the immune system responds to vaccination (Hazan et al., 2023;Cermakian et al., 2021;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.
The molecular basis of the immune-circadian interaction has been previously linked to specific processes within the innate immune system (Labrecque and Cermakian, 2015;Jerigova et al., 2022).Indeed, molecular clocks have been characterized in macrophages, neutrophils, eosinophils and natural killer (NK) cells (Keller et al., 2009).For example, 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 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 Per2 in lymph node (Fortier et al., 2011;Nobis et al., 2019a) and circulating naïve CD8 + T cells (Dimitrov et al., 2009), the molecular source(s) of this daily variation remains unknown.Additionally, this circadian variation in CD8 + T cell activity may, in part, depend on cell-autonomous mechanisms, such as clocks within those cells (Nobis et al., 2019a;Wang et al., 2023).Through transcriptomic analysis, a recent study identified several genes within the CD8 + T cell activation pathway that exhibited rhythmicity in their expression in naïve CD8 + T cells, and subsequently demonstrated that certain molecular players within this activation pathway, such as ZAP70, AKT and mTOR, demonstrated an increased activation following a day-time vaccination, as well as an increased ability to combat an infectious challenge after a day-time vaccination (both compared to night-time vaccination) (Nobis et al., 2019b).Moreover, another study showed an increased level of ZAP70 prior to vaccination in the day-time (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 a vaccine input and its underlying dynamics.We adapt 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 TCR phosphorylation.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 activation postvaccination 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 on a selected parameter to develop a deeper understanding of the dynamics of this modulation.We then expand this approach to a multi-parameter 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 T cell gene regulatory 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 this circadian modulation, and that the latter is dependent on both the number of phosphorylated T cell receptors (TCRs), which we term the vaccine load, presented to the cell, 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 the mTORC1 protein complex, may be independently regulated by circadian rhythms, pointing to the existence of a complex, multi-layered circadian influence on CD8 + T cell signaling.Our population modeling further expand on these ideas using a multi-parameter framework to produce more physio-logical conclusions about the influence of a circadian rhythm on CD8 + T cell activation.

Gene regulatory network model
To explore the dynamics of circadian influence 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 + T 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 a 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).The rate equations and parameter values of this model are provided in Table 6.1.

The three key axes of the model
To explore the impact of a 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 green 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 a 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 axis is more sensitive to circadian-induced perturbations, helping us identify potential molecular targets of this circadian rhythm.

Periodic forcing by circadian rhythms
To investigate the effect(s) of a 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 variation 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 a 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 α, representing 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 circadian rhythms is consistent with how it is routinely characterized in experimental literature (Elliott et al., 1972;Hickey et al., 1984;Naitoh et al., 1985).
To account for physiological limitations of how much a parameter can vary (see Appendix 6.3.2),we assume that α p ∈ [0, 0.5].The periodic forcing defined in Eqn. 1 produces oscillations with a 24 h period and a phase shift of 0 for a positive feedback parameter and π for a negative feedback parameter; this reproduces the effect of a circadian rhythm, and leads to a variation in the signaling molecule associated with p t , depending on time of day.Increasing the circadian coefficient α causes the amplitude of p t to increase and thus the range of values of p t to expand in a parameter-specific manner.
3. Outcomes of the single cell model

Excluding the oscillatory regime
The original parameterization for the ZAP axis described by Zheng et al. (2005)  Identifying the parameters responsible for producing the key features of model dynamics will provide us with potential targets for producing the day-night variation in model output.By performing a dynamical systems analysis of the model, one can predict how parameter perturbations, due to the CD8 + T cell circadian rhythm, can modify model behaviour.To accomplish this, we make the assumption here that a circadian rhythm on CD8 + T cell signaling would primarily affect the rate activation kinetics of signaling molecules within the CD8 + T cell activation pathway.This is supported by the day-night variation in active forms of CD8 + T cell signaling molecules shown in Nobis et al. (Nobis et al., 2019b).
We begin first by performing parameter sweeps and bifurcation analyses within each axis of the gene regulatory network model to identify the broad signaling pathways defining the dynamics of our model.Once a 'key' axis is found, individual parameters within the axis are further explored using both parameter sensitivity and threshold analyses to identify targets for this circadian rhythm (i.e., periodic forcing).The stable and unstable branches are linked by Hopf bifurcations (HB1 and HB2) where envelopes of unstable limit cycles (blue dashed lines) emerge and terminate at homoclinic bifurcations (HM1 and HM2, respectively).As k10f decreases, a jump in the ZAP active steady state occurs upon crossing the HB2 bifurcation point.We also note that there is a bistability region within which the two branches of stable steady states overlap (shown by the region where the two red solid lines cover the same range of T CRp parameter values).This region is bounded by HB1 and HB2.Since the Hopf bifurcations produce envelopes of unstable limit cycles, periodic solutions are not observed when transitioning from one stable branch of steady states to the other.Additionally, since the bistable regime is bounded by the two Hopf bifurcations HB1 and HB2, rather than the two saddle nodes SN1 and SN2, as typically seen in traditional bistable switches, the bistable regime and the hysteresis of the model are quite limited in range.This bifurcation diagram suggests that certain molecular players within the ZAP axis are quite important in modulating its response to a vaccine input; however, it remains unclear which specific parameters within the axis are most sensitive to perturbations in their values.
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 These results thus suggest 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.

Parameters of ZAP 70 and SHP 1 activation are the most sensitive to perturbation
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 i within the ZAP axis, and is done using a MATLAB package written by Marino et al. (Marino et al., 2008), which incorporates Latin Hypercube Sampling (LHS) 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., 2019b), 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 is the primary negative feed- back parameter implicated by this analysis, it will be our target parameter in the single-parameter analysis conducted in the following section.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, a 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 this circadian rhythm can cause periodic oscillations in ZAP active with large amplitudes.This then poses the question of how such mechanism interacts with vaccines?
To answer this question, we will investigate how the impact of a circadian rhythm on CD8 + T cell activation is affected by the T CRp (the number of phosphorylated TCRs), which we assume to be an indicator of vaccine strength.Thus, an elevated T CRp is assumed to demonstrate a vaccine presentation to the CD8 + T cell, and a higher vaccine load is similarly assumed to trigger a higher number of phosphorylated TCRs, thus leading to a larger T CRp model input.This is first done by considering two other values of T CRp: one that is higher than the one used in Fig. 4 Taken together, these results demonstrate that the day-night variation 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.

Impact of circadian coupling on temporal evolution of CD8 + T cell signaling
Now that we have determined the underlying steady state dynamics of circadian modulation of CD8 + T cell activation with different vaccine loads, we 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 (pre-vaccination) with an initial amount of T CRp given by T CRp 0 = 100 mol; a nonzero basal level of T RCp 0 is chosen for this state S 0 to account for 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 ), i.e., when the CD8 + T cell has been presented with a vaccine.The effect of a 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 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 a vaccine is applied at both daytime (at t = 12 hr) and nighttime

Characterizing the filtering properties of the ZAP axis
Throughout this study, we have explored the implications of a 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 a circadian rhythm acting directly on the ZAP axis is transmitted to the downstream LAT and CD28 axes within the model in the absence of a vaccine input, i.e., in the S 0 state.Doing so would indicate whether a circadian rhythm within the ZAP axis is capable of producing a day-night variation in the signaling of other downstream molecular players associated with these two axes.If not, this may suggest that those molecular players may be independently regulated by the circadian clock.
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. 1 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.Since the model is in S 0 , this value of α k10f is chosen as it creates the largest variation in steady state amplitude, allowing us to more easily visualize how a circadian rhythm is filtered through the model, without the possibility of jumping to an upper, 'activated', steady state branch.
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 for 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., 2019b), despite producing them in the other signaling molecules of this axis.
This was resolved by imposing a periodic forcing component, representing the influence of a 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 variation (Nobis et al., 2019b), these modeling results suggest that a circadian rhythm within the ZAP axis may not be sufficient to explain these downstream oscillations, pointing to mT ORC1 and its related players being secondary, independent targets of the molecular circadian clock.

Methodology
To better understand CD8 + T cell responses to circadian oscillations and vaccines, it is imperative to obtain a broader view of model dynamics at the T cell population level.We do so by simulating the steady state dynamics of n = 1000 model realizations of the ZAP axis, each with a randomly selected parameterization of the model.Similarly to the single-cell approach, the oscillatory regime of the model was excluded using the framework described in Appendix 6.3.This has allowed us to confine parameter sampling of each realization to the non-oscillatory regime only.

One-parameter macro-bifurcation analysis
We begin our population level analysis by first targeting one specific parameter, namely, k10f k .This is done by randomly assigning values for this parameter within the range [0.5 × k10f 0 , 1.5 × k10f 0 ] for each model realization k, where k10f 0 denotes the default parameter value of k10f listed in Table 6.ZAP active is purely elevated, regardless of the value of k10f t .

Multi-parameter macro-bifurcation analysis
Unique to the macro-bifurcation analysis of the previous section is the ability to expand it to examining n parameters simultaneously, an advantage not afforded by the standard bifurcation analysis that is typically limited to at most 2 parameters (as in Fig. 4E).To harness this advantage, we now include variability not only in k10f , but also in the positive feedback parameters k4f, k9f 1, negative feedback parameters k7f, k11f and our manually derived threshold parameter k8f 1 (see Appendix 6.3.2).We have chosen these parameters because they can act as possible molecular targets of a 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 and k8f 1, please see Tables 6.1).Combining their effects in a macro-bifurcation diagram can thus provide a more physiological understanding of how a circadian rhythm can impact CD8 + T cell activation.Simulating 1000 realizations of the model, by randomly, and independently, 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 Further analysis, coupled with experimental work, may help determine the mechanism by which these additional rate parameters interact with the circadian rhythm directly to produce this increased sensitivity, providing a more complete understanding of the complex population-level immune-circadian interaction.

Discussion
In this study, we provided a detailed analysis of the dynamics of the signaling molecules involved in producing the day-night variation seen in CD8 + T cells signaling.Additionally, we explored how this day-night variation in CD8+ T cells downstream of the T cell receptor (TCR) may be altered by different levels of TCR phosphorylation.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 a circadian rhythm within the ZAP axis, a global sensitivity analysis was conducted using Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) (Marino et al., 2008).This was done to identify how parameter perturbations within the ZAP axis can affect the dynamics of its output signal, namely ZAP active .Since this 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 and k8f 2 were the most sensitive to perturbations, conclusions that are aligned with the results of transcriptomic analyses of CD8 + T cells in mice (Nobis et al., 2019b).
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 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 variation.This day-night variation was also shown to be affected by the quantity of TCR phosphorylation, which we designate as the vaccine load, via the input parameter T CRp.This suggests a dependency between the level of TCR phosphorylation after vaccine input and the influence of this circadian rhythm, in such a way that there exists a threshold in T CRp below of which a circadian rhythm may not affect CD8 + T cell signaling.These results thus show that there is a potential complex relationship between the vaccine load presented to a CD8 + T cell and the molecular circadian clock.
Through a time-series analysis at 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 variation and that such a circadian effect ceases when a vaccine is applied.This suggests that the physiological circadian rhythm acts during the baseline, 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., 2019b).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 revealed 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 and 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 variation was experimentally observed in mT ORC1 * (Nobis et al., 2019b), our results suggest that this day-night variation may be due to independent circadian rhythms acting directly on mT ORC1, rather than being due to the activity of an upstream circadian rhythm within the ZAP axis.
To produce a more physiological representation of CD8 + T cell responses to vaccination and how that is impacted by a circadian rhythm, we implemented a multi-dimensional population level modeling approach.This allowed us to explore the effect of a CD8 + T cell circadian rhythm on multiple molecular players simultaneously, through what we termed as a macro-bifurcation diagram.In these simulations, key circadian-dependent parameters were randomly sampled within their respective ranges, and produced three regions of activity with respect to vaccine load, defined by T CRp: insensitive, sensitive (macro-bistable) and elevated.Additionally, we demonstrated that adding additional rate activation parameters to the list of circadian-dependent parameters expanded the sensitive region of the macro-bifurcation diagram.This suggests that the increased variability in the number of circadian-dependent parameters increases the sensitivity of a population of CD8 + T cells to vaccination (by lowering the minimum T CRp needed for CD8 + T cells to be in the elevated steady state).That said, further analysis should be done to link this increased sensitivity to the circadian rhythm directly, by identifying which parameter sets bias T cells to the upper steady states, providing us a more complete understanding of the specifics of this interaction on a population level.
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 a circadian rhythm on CD8 + T cell signaling, in the context of vaccination, and serves a valid framework for assessing the influence of circadian rhythms in other models.

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

Variable Differential Equation
ZAP axis CD28 axis

Rate Equation Components
ZAP axis      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, k10f (which represents SHP1 activation) and k8f 1 (which represents ZAP activation); the latter is 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, chosen from within this range.Since we are interested in the effects of a CD8 + T cell circadian rhythm, this parameter range should include the degree of modulation due to this 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 n 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 this 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., 2019b).We will therefore use this range to describe the influence of a circadian rhythm on CD8 + T cell signaling.To produce this p range, we choose the range of a viable circadian coefficient α to be α ∈ [0, 0.5].It should be noted that this only defines the maximum range of p; some realizations, k = 1, 2, . . ., n, 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 ] variation in cytokine production can be traced linearly to a variation 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 variation 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 dynamic behaviour
With n model realizations selected from parameters within physiological ranges, each individual realization, k = 1, 2, . . ., n, can be integrated and the steady-state solution associated with each one 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 1) 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 n model realizations, each with its unique parameterization of the m circadian-affected parameters, and with a label designating the classification of each solution k = 1, 2, . . ., n 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.Each model realization is linked to its respective steady state dynamics according to the values of its parameters of interest.Such visualization 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.

Figure 2 :
Figure 2: (A) Time series simulations of ZAP active , the output of the ZAP axis, at different values of k10f and k4f specified in the legend, where k10f and k4f are the parameters for SHP 1 phosphorylation, SF K * and ZAP p binding, respectively.(B) One-parameter bifurcation of ZAP active with respect to k10f .Red (gray) solid (dashed) lines represent stable (unstable) branches of equilibria, while blue dashed lines represent envelopes of unstable limit cycles emerging from the Hopf bifurcations: HB1 and HB2, and terminating at the homoclinic bifurcations: HM1 and HM2.SN1 and SN2 are saddle-node bifurcations.(C) Time series simulations of P LCγp as a percentage relative to its total concentration: P LCγ + P LCγp, at different values of k12f and k13f specified in the legend, where k12f and k13f are the forward rate constants of LAT and P LCγ phosphorylation, respectively.(D) One-parameter bifurcation of P LCγp percentage with respect to k12f .Red line represents branch of stable equilibria.(E) Time series simulations of mT ORC1 * as a percentage relative to its total concentration: mT ORC1 + mT ORC1 * , at different values of k42f and k43f specified in the legend, where k42f and k43f are the forward rate constants of RhebGDP and mT ORC1 activation, respectively.(D) One-parameter bifurcation of mT ORC1 * percentage with respect to k42f .Red line represents branch of stable equilibria.All time series simulations and bifurcation diagrams are plotted with T CRp = 27000 mol; for all other parameter values, seeTable 6.1.
and mT ORC1 * percentages (relative to their individual total concentrations) 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, respectively) 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 * percentages 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.

Figure 3 :
Figure 3: 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 and 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|.

3. 3 .
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)

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) One-parameter bifurcation of ZAP active with respect to k10f at 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 extrema identified in (A) are superimposed as yellow and light blue dots.
(A), namely, T CRp = 40000 mol, and another that is lower, i.e., T CRp = 10000 mol, and then plotting the one-parameter bifurcations of ZAP active with respect to k10f at these two T RCp values as shown in Fig.4(C,D), respectively.By superimposing the orange and blue dots representing the day and night extrema of k10f t when α k10f = 0.5 on these two bifurcation diagrams, we see that in both realizations, no Hopf or saddle node bifurcations end up separating the two dots to allow for drastic time-of-day dependent changes in ZAP active seen in 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 this circadian rhythm can produce the largest day-night variation 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) demarcates 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 of T CRp is varied, indicating that there is an intermediate range of T CRp within which the largest day-night variation in ZAP active can be generated.This means that the level of TCR phosphorylation, or level of vaccine load, can alter T CRp and cause the regime of bistability to shift accordingly.Now that we have explored how differing vaccine loads may affect bistability with respect to T CRp and k10f , we now examine how changing the strength of this circadian rhythm through α k10f can affect outcomes.By tak-ing 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.According to this model 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 a realization thus indicates that if the influence of this 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 .
(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 TCRs become phosphorylated, due to vaccine application, 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.

(
at t = 24 hr) as shown inFig.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 back to the lower branch.These results thus corroborate 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 a circadian rhythm on early events of CD8 + T cell signaling is significantly attenuated post vaccination.Based on these observations, we conclude that this circadian rhythm primes or facilitates the activation of naïve CD8 + T cells to produce the day-night variation in ZAP active output, rather than directly controlling it.

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 pre and post vaccination applied at nighttime (t = 24 hr).Inset: magnification of the time series simulations prio to vaccination to highlight the small amplitude oscillations associated with the state S 0 .(C)Time series simulations of ZAP active pre 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 pre and post vaccination applied at nighttime (t = 24 hr).(F) Time series simulations of ZAP active pre 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, (B) LAT, and (C) CD28 axes when k10f t is oscillatory with α k10f = 0.5.The larger the ω i value, the stronger the influence of the circadianinduced 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.
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 this graph represents a steady state value of ZAP active at a given k10f k within the previously defined range.Doing so provides a snapshot of 1000 steady states of the model, each at a different, randomly selected, circadian time defined by the value of k10f k ∈ [0.5 × k10f 0 , 1.5 × k10f 0 ] (recall that we assumed that this circadian rhythm affects the value of k10f when we defined k10f t ).Plotting the steady states of these 1000 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 this circadian rhythm on T cell activation through k10f , but also the effects of vaccine strength through T RCp.We refer to this plot as a 'macro-bifurcation' diagram.
), and plotting the steady state values of ZAP active produced by the model with respect to T CRp produces the macro-bifurcation diagram shown in Fig. 7(B).Although the distinction between insensitive, sensitive and elevated regions remains regardless of the number of parameters included, the macro-bistable region of the model (the sensitive region) becomes more expanded (i.e., more sensitive) compared to that shown in Fig. 7(A): T CRp ∈ [2, 3] • 10 4 mol in panel (A) vs T CRp ∈ [1.5, 3.5] • 10 4 in panel (B).This demonstrates that our theoretical framework is maintained when adding additional parameters, and suggests that incorporating variability in additional parameters in the macro-bifurcation analysis increases the range of the macro-bistability by introducing more heterogeneity in CD8 + T cells.This increased heterogeneity due to the circadian rhythm may improve the sensitivity of CD8 + T cells at the population level, allowing them to activate at a lower vaccine load.
ZAP70 and mTOR translate to changes in CD8 + T cell activation.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 .
Oscillatory dynamics of the model 895 With the parameterization of the model described in Zheng et al. (2005) 896 and Perley et al. (2014), the model is able to produce intrinsic oscillatory 897 signals, shown in Fig. S1(A).For instance, setting T CRp = 20000 mol 898 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(B), produces a diagram that effectively eliminates thebistable behavior of the model and produces an oscillatory regime defined by the envelopes of stable limit cycles that originate from 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 with periodic spiking in ZAP active similar to that shown in Fig.S1(A).This spiking pattern in ZAP active has not been observed experimentally, indicating that this is a non-viable parameter regime of the model.Physiologically, 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 aninvalid 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.6.3.Eliminating the intrinsic oscillatory regimeThis approach is adapted from the Dynamics Visualisation based Parallel Coordinates (DYVIPAC) framework byNguyen et al. (2015).Broadly, it involves a population-based approach to uncover system's dynamics.The m parameters whose effects on model dynamics are the main focus of the analysis are first selected, and a range of values for each parameter is set.A sample of n model realizations (i.e.simulated steady state responses) is then created, each with a randomly selected parameterization from within predetermined ranges.These realizations are integrated starting from the same initial conditions, and their steady state solutions are classified according to outcomes.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

[ 2 .
72, 6.12] (mol • s) −1 .Applying parallel coordinate plot on this new regime now only produces two outcomes: the inactive and elevated states, as shown in FigureS2(C).Through this new k8f 1 parameterization, we observe that the oscillatory regime is indeed eliminated, while still preserving both the inactive and elevated model solution.Thus, this new model parameterization can safely allow us to simulate the effects of a circadian rhythm on CD8 + T cell signaling impacting directly k10f and k8f 1 without worrying about intrinsic oscillatory dynamics that could be produced by the model; this is because we are certain that any combination of these three parameters under this parameterization considered in FigureS2(C) would not produce the non-physiological oscillatory solutions.This constrained parameter regime of k8f 1 is used in both our single-cell and population-level modeling simulations to avoid the emergence of oscillatory solutions.

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 yaxis.Each segmented colored line represents a given model realization (i.e., steady state response) sampled randomly from within a pre-determined range.The segmented lines are color-coded 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).Range of parameter values for (A) & (B): T CRp ∈ [0, 30000], k8f 1 ∈ [4.17, 12.50], k10f ∈ [4.32, 12.95].(C) The analysis repeated with a constrained parameter range for k8f 1.Notice the lack of yellow oscillatory solutions.Range of parameter values for (C): T CRp ∈ [0, 30000], k8f 1 ∈ [2, 6], k10f ∈ [4.32, 12.95].

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

Table 2 :
Perley et al. (2014)f the various signaling molecules within the model.Equations in the 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).