Extracting Circadian and Sleep Parameters from Longitudinal Data in Schizophrenia for the Design of Pragmatic Light Interventions

Abstract Sleep and circadian rhythm dysfunction is prevalent in schizophrenia, is associated with distress and poorer clinical status, yet remains an under-recognized therapeutic target. The development of new therapies requires the identification of the primary drivers of these abnormalities. Understanding of the regulation of sleep–wake timing is now sufficiently advanced for mathematical model-based analyses to identify the relative contribution of endogenous circadian processes, behavioral or environmental influences on sleep-wake disturbance and guide the development of personalized treatments. Here, we have elucidated factors underlying disturbed sleep-wake timing by applying a predictive mathematical model for the interaction of light and the circadian and homeostatic regulation of sleep to actigraphy, light, and melatonin profiles from 20 schizophrenia patients and 21 age-matched healthy unemployed controls, and designed interventions which restored sleep-circadian function. Compared to controls, those with schizophrenia slept longer, had more variable sleep timing, and received significantly fewer hours of bright light (light > 500 lux), which was associated with greater variance in sleep timing. Combining the model with the objective data revealed that non 24-h sleep could be best explained by reduced light exposure rather than differences in intrinsic circadian period. Modeling implied that late sleep offset and non 24-h sleep timing in schizophrenia can be normalized by changes in environmental light–dark profiles, without imposing major lifestyle changes. Aberrant timing and intensity of light exposure patterns are likely causal factors in sleep timing disturbances in schizophrenia. Implementing our new model-data framework in clinical practice could deliver personalized and acceptable light–dark interventions that normalize sleep-wake timing.


Further participant details Living situations
At the time of the study, thirteen community-living patients lived on their own and 15 of the controls; one patient and 1 control each lived with a partner; one patient woman lived with a child and 1 control woman with two children; five single patients lived in shared accommodation (2 with parents, 2 in supported housing, 1 in hospital awaiting accommodation) and 4 single controls lived in shared flats (2-and-2).

Clinical status and medication
Activities of daily living and social activity by self-report and notes from mental health services was rated using the General Assessment of Function scale (GAF, Frances et al, 1994) and their score ranged from 34 to 70 (median 50) [1]. Behaviourally controlled, persistent hallucination or delusions were present in eight individuals. Regularly prescribed psychotropic medication to manage behavioural disturbances included antipsychotics only (n = 12), additional psychotropic medication to control extrapyramidal (antimuscarinics, n = 3) and seizure-inducing side effects of antipsychotics (anticonvulsants, n = 4). Depressive or anxiety symptoms were experienced by eight people at some point during their illness, two of these were treated with a selective serotonin inhibitor. Antipsychotic medication and dose (ranked by chlorpromazine equivalent), age, gender, length of illness, drug delivery mode and absence/presence of positive symptoms are detailed in Table S1.
For those with schizophrenia, illness duration ranged from 2 to 33 years (median 10 years).
The patients were not high functioning but were average individuals with histories comparable to [2].
In [1], the presence of positive symptoms, age, length of illness, GAF and antipsychotic drug dose were not found to correlate with measures of sleep, rest-activity and melatonin, but differences in melatonin levels between individuals with and without positive symptoms were an effect of age. In our current study, we additionally found no correlation between measures of light and chlorpromazine equivalent (Spearman's rank: hours of bright light, p = 0.97; total light, p = 0.74; total log light p = 0.50).

Monitoring of sleep-wake timing
In view of the functional disability associated with schizophrenia, to enhance data quality weekly home visits were made by the investigator. At these visits, actigraphy, light and diary entries were compared and any discrepancies discussed with the participant. This high level of supervision may not be practical in a clinical setting. However, in another recent study, we found surprisingly high rates of adherence with remote rest-activity monitoring in populations with chronic schizophrenia, some with high symptom burden of positive and negative symptoms [3]. We therefore believe that our approach is generalisable.

Further statistical details and results
Linear mixed effect models were used to investigate the repeated measures of sleep duration, sleep onset and offset times. In all cases, fixed effects were participant group (control or schizophrenia) and gender. Participant was included as a random effect.
Results are summarised in Table S2. The effect of group on sleep measures is reported in the main text. We found no significant effect of gender.

Mathematical modelling
The model describes sleep and wake states as the firing of sleep-promoting and wake-promoting populations of neurons respectively. Switching between sleep and wake occurs as a result of homeostatic and circadian drives to sleep-promoting neurons. The homeostatic sleep drive increases with time awake, as the physiological pressure for sleep increases. The circadian drive models oscillatory input from the master circadian clock.
During the biological day the circadian clock promotes wakefulness whereas over the biological night it promotes sleep [4][5][6].
In addition, the model captures the impact of light exposure and endogenous circadian

Equations
As in [4,5], the interaction between sleep promoting and wake promoting neurons was described by equations for their mean electric potential, V v and V m , respectively, Here, describes the sigmoidal relationship between the potential and the firing rate of the neurons, Q v,m . The parameter Q max is the maximum possible firing rate; θ is the value of the potential V v,m at Q v,m = Q max /2 and σ determines the width of the sigmoid.
The parameter τ gives the typical timescale of the neuronal process and the parameters ν vm,mv weight the input from population m to v and v to m respectively.
The D v,m are 'drives'. The drive to the wake promoting neurons, D m was fixed. The drive to the sleep promoting neurons, D v , consisted of homeostatic and circadian components, where H(t) describes the homeostatic sleep pressure and C(t) represents the circadian wake promoting rhythm. As discussed further below, C(t) is approximately sinusoidal with amplitude close to one, so ν vc gives the amplitude of the circadian wake propensity rhythm.
The term A v is a background excitatory input to the sleep promoting neurons.
The homeostatic process is given by Here, H/χ is the rate of removal of somnogenic chemical and µ Q m /χ is the rate of production.
The circadian rhythm C(t) was modelled by a forced van der Pol oscillator [8], where the forcing represents the light intensity dependent signal to the suprachiasmatic nucleus from photoreceptors in the eye. The model was originally constructed to accurately replicate human phase response data [8]. Specifically, where Here, x and y are the variables for the van der Pol oscillator and n is the fraction of activated photoreceptors. The term B describes the forcing that occurs as a result of the action of light of intensityĨ on photoreceptors. The magnitude of B is dependent on phase, as modelled by the term (1 − bx)(1 − by) in equation (9). The rate of saturation of photoreceptors is given by λα 0 Ĩ I 0 p and the rate for decay is given by λβ, where the values of λ, α 0 , I 0 , p and β have been determined from experimental data [7]. The parameter τ c is the intrinsic period, γ is the stiffness of the oscillator.
Note that the intrinsic period of the van der Pol oscillator was dependent on light intensity since for diurnal species intrinsic period is shorter for higher light intensities (Aschoff's rule). This dependence of the intrinsic period on light intensity is controlled by the parameter k in equation (6).
When light intensity I(t) is provided by measurements in the field, we make the assumption that the relationship between light intensityĨ incident on photoreceptors in the eye and measured light isĨ = I(t). In this case, information on sleep timing is implicitly given within the light data since lights tend to be turned off during sleep.
When using the model in predictive mode, we make no assumption that light is turned off at night and light is 'available' at all times of day. In these cases,Ĩ is modelled as where H is the Heaviside function. This means that light impinging on photoreceptors is zero if the firing rate of wake promoting neurons, Q m , is greater than a threshold value Q th i.e. light is turned off during model-predicted sleep. See Fig. S1(a).
As in [5], the circadian wake propensity rhythm C(t), was modelled as   Table S3. The sleep/wake regulation and circadian parameters are the same as those used in the original version of the model [4].
Equations were integrated using the stiff solver ODE15s in MATLAB [9]. A typical solution for equations (1)-(11) is shown in Fig. S1.  Table S3. This figure is reproduced from the Supplementary Material for [5].

Modelling changes in sleep duration and timing
In equations (1) -(5), sleep duration is determined by three separate groups of parameters.
First, parameters which relate to the rate of increase of homeostatic sleep pressure (µ and χ). Second, by parameters which relate to the thresholds at which switching between wake and sleep occur (A v and A m ) and finally, by ν vh which determines the sensitivity to homeostatic sleep pressure. We note that there is a scaling between ν vh and µ which means that mathematically these two parameters are equivalent: one of them could be scaled from the equations. Here, as in previous publications we retain both.
Previously, we have modelled changes in sleep duration across the lifespan by changing the parameter µ [10]. However, studies of slow wave sleep provide no evidence that those with schizophrenia have a higher sleep need [11]. Hence here we considered the parameters that govern wake propensity, namely A v and A m . Wake propensity could be reduced by lowering the baseline input to wake promoting neurons (A v ). However, paradoxically, with the form we have taken for the rate of production of somnogenic chemicals in equation (5), reducing the drive for wake results in a decrease in sleep duration. This is because of a 'use dependent' effect: reducing the drive for wake also reduces the firing rate of the wake promoting neurons and thus the build up of homeostatic sleep pressure during wake. Hence, here we have modelled reduced wake propensity by reducing the mean level of the input to sleep promoting neurons, A v in equation (4). This is equivalent to lowering both thresholds of the wake propensity rhythm in the two process model [12][13][14] and also Fig. S2.
Parameters that control sleep timing include the intrinsic circadian period τ c in equation (7) and the circadian amplitude, ν vc . Here, since there are known individual differences in intrinsic circadian period, [15], we elected to vary intrinsic circadian period, see Fig. S2.

Parameter fitting
For each individual, we selected the value of A v that matched mean sleep duration in the model to mean observed sleep duration. Since there is a monotonic relation between sleep duration and A v , see Fig. S2, this is straightforward.
Next, for each individual, we fed the raw light data into the mathematical model and  Bottom right, wake propensity as measured by deviation from the default value for A v , with negative numbers indicating decreased wake propensity i.e. wake propensity = A v + 10.2. Note, here for consistency with the upper panels, the horizontal axis is scaled so that lower values of wake propensity, hence longer sleep duration, is to the right.
Performing the fitting sequentially, as described, is possible because changing intrinsic period has minimal effect on sleep duration. It is computationally more efficient that simultaneously fitting for both parameters.
We investigated whether there were gender differences in the distributions of fitted parameters by using the Wilcoxon rank-sum test. For fitted wake propensity, there was no significant difference between men and women (p = 0.409). In line with previous observations [15], fitted intrinsic circadian period was on average slightly longer in men than in women, but did not reach significance in our sample (p = 0.051).

Modelling light interventions
Motivated by the shape of real world light profiles [16,17] and as introduced in our previous work [5], we used the function This ramps up to a maximum level l i 1 during core day light hours on day i and reduces to a baseline of l 2 during the evening, where the transition from l 2 to l 1 occurs around t = s 1 and from l 1 to l 2 around t = s 2 . The speed of the transition is determined by the parameter c. For appropriate parameter choices, the profile given by equation (12) gives a good match to the mean hourly light intensity profile observed in the winter by participant sz01 who exhibited non-24 h rhythms, see Fig. S3.
In order to systematically investigate the effect of different light profiles, the values for the 'evening light' l 2 and the 'day light' l 2 were varied over a range of relevant values. Four different combinations of the physiological parameters were used to demonstrate the impact of intrinsic period τ c and wake propensity A v on the regions of entrainment. The results are shown in Fig. S4. As one might expect, these show that increasing intrinsic circadian period makes it harder to entrain to the 24 h day. Thus those with a long intrinsic period would need higher levels of light during the day and/or lower levels of light in the evening to achieve a desired spontaneous wake time than those with an intrinsic circadian period closer to 24 h. Decreased wake propensity can either make it easier or harder, depending on the intrinsic period. Note that the longer sleep duration means that a reduced wake propensity makes it harder to get up early in the morning.
The interaction between wake propensity and intrinsic period means that without the mathematical model it is hard to predict the outcome of a change in available light.
We note that the overall shape of the entrainment regions is consistent with the way that the Arnold entrainment tongues [5] shift as a function of evening light and day light.