The Neurospora photoreceptor VIVID exerts negative and positive control on light sensing to achieve adaptation

The light response in Neurospora is mediated by the photoreceptor and circadian transcription factor White Collar Complex (WCC). The expression rate of the WCC target genes adapts in daylight and remains refractory to moonlight, despite the extraordinary light sensitivity of the WCC. To explain this photoadaptation, feedback inhibition by the WCC interaction partner VIVID (VVD) has been invoked. Here we show through data-driven mathematical modeling that VVD allows Neurospora to detect relative changes in light intensity. To achieve this behavior, VVD acts as an inhibitor of WCC-driven gene expression and, at the same time, as a positive regulator that maintains the responsiveness of the photosystem. Our data indicate that this paradoxical function is realized by a futile cycle that involves the light-induced sequestration of active WCC by VVD and the replenishment of the activatable WCC pool through the decay of the photoactivated state. Our quantitative study uncovers a novel network motif for achieving sensory adaptation and defines a core input module of the circadian clock in Neurospora.


Table of Contents
. Dilution series to test the probe efficiencies. Figure S2. VVD stabilizes WCC protein levels Figure S3. Transfer from light to dark results in the onset of oscillations. Text S1. Mathematical Model Text S2. Futile cycling provides the means for repeated induction References MATLB code for all model simulations are available as separate supplementary files.

Supplementary Figures
Figure S1  Figure S1: Dilution series to test the probe efficiencies.
To ascertain the relative differences between vvd, wc-1 and frq RNA levels and thus be able to compare them with each other, we tested the efficiency of their respective primers and probes in a qRT-PCR using a dilution series of wildtype gDNA as template. Calculating the delta CT with vvd at the lowest DNA template concentration as reference, we determined the fold difference in the probe efficiencies: vvd = 1; frq = 2.086 and wc-1 = 1.552. The values of the experiment (Fig. 1) were then corrected accordingly, with reference to vvd. qRT-PCR with a dilution series (6.25, 25 and 100 ng) of wild type gDNA as template probing for frq (diamonds), wc-1 (triangles) and vvd     Self-sustained oscillations in constant darkness are thought to arise through delayed negative feedback of the WCC target FRQ on the activity of the WCC (Brunner and Káldi, 2008;Baker et al , 2012). These delays are included in the form of a series of intermediate steps between FRQ synthesis and FRQ-dependent inhibition of the WCC: − − → ∅, where the final FRQ mediates phosphorylation of the WCC. The additional parameters are: kp1 = 0.375 hr −1 ; kp2 = 0.375 hr −1 ; kpn = 0.375 hr −1 A small set of modifications had to be made to the other parameters to ensure that adaptation to light is retained. Modified parameter values: k3=10.0 hr −1 ; K3 = 6.0 au; n3 = 5.0; K4 = 61.0 au; k−4 = 0.3 hr −1 ; k11 = 43.0 hr −1 ; K11 = 0.07 au; k19 = 0.24 hr −1 ; q6 = 750.0 hr −1 ; Q6 = 0.07 au. The modification of parameter values and the introduction of time-delayed FRQ kinetics do not affect the adaptation behavior of the system. In all cases, the blue curves correspond to the wild type, the red curves to the vvd-mutant. Solid curves represent the model simulation results; dashed lines, the experimental results. (A) Simulation of light induction of frq mRNA followed by transfer to darkness shows the characteristic adaptation to light before the onset of oscillations in the dark. In the vvd-, frq mRNA levels decrease slower than in the wild type after the transfer from light to dark (cf Figure 6) resulting in a delay in onset of oscillations. Apart from this, the period and amplitude of the oscillations are the same between the wild type and the vvd-mutant, confirming that oscillations in free-running conditions are independent of VVD. (B) frq mRNA levels oscillate in the dark with oscillations phase-shifted in the vvd-by approximately 4 hours, an observation that has previously been noted (Heintzen et al , 2001;Elvin et al , 2005;Hunt et al , 2007). (C)-(D) Responsiveness and adaptation to increasing light steps show that the wild type remains responsive and adapts within 4 hours, while the vvd-mutant requires longer for adaptation to be achieved. (E)-(F) A light to dark transfer shows that transcriptional activity for vvd and frq mRNA are unregulated for up to 4 hours in the vvd-mutant, while activity drops within 1 hour in the wild type. (G)-(H) Steady-state vvd mRNA levels are shown experimentally to increase with increasing light intensity, while the levels remain constant in the vvd-mutant. Model simulations agree with the experimental results.
Supplementary Text 1

Mathematical model
Based on time-resolved experimental data we constructed a mathematical model to understand the role of VVD as both a negative and positive regulator in light adaptation.
We outline below, some of the model features describing the interaction between the circadian transcription factor the White Collar Complex (WCC) and its negative regulators VIVD (VVD) and Frequency (FRQ). Saturation terms have been used for enzymatic reactions: induction of gene expression and phosphorylation steps by FRQ. Photoactivation, complex formation, dissociation and (non-enzymatic) photoadduct have been modeled by mass action terms.

Light-activation of WCC and VVD
WC-1 and VVD contain a flavin-binding light-oxygen-voltage (LOV) domain, the blue-light sensors in the light receptors of plants and fungi. Light triggers a LOV-mediated homodimerization of the WCC, which induces expression of a large number of genes.
We model the light activation of both the WCC and VVD by a two-step process describing first activation by light (to give the light-activated form of the WCC, WCC*) followed by homodimerizaton of WCC*. Both steps are modeled using mass action kinetics.
The effect of light is modeled through a direct change to the activation rates of the dark forms of unphosphorylated and phosphorylated WCC (WCC, WCC p ), k 1 , and VVD (VVD), q 1 , to form the corresponding light-activated forms, WCC*, WCC* p and VVD*.
Values for the activation rate are estimated from Malzahn et al , 2010 and given in Table S1.

Light intensity
Value Table S1: Light-activation rates, k1, q1, for each of the experimental light intensities used.

Transcription
The WCC drives transcription of many genes, including vvd, frq and wc-1.
In the dark, the inactive, dark form of the WCC can drive expression of frq by binding to the clockbox (C-box) element of the frq promoter (Froehlich et al , 2003). This is given by the second term in the equation for d [FRQ] dt , where WCC binds with affinity K 14 . Additionally, light-activated WCC (in the form of the homodimer, WCC*-WCC*) can also bind to the C-box and to the proximal lightresponsive element (LRE) of frq. Binding affinity of WCC*-WCC* to the C-box is given by K 14 ′ , while WCC*-WCC* binds to the LRE with affinity K 13 .
vvd transcription is driven only by the light-activated homodimer (WCC*-WCC*) with a maximal transcription rate of q 6 .
Transcription of wc-1 in the dark is independent of the WCC, given by the basal rate k 12 . In light, transcription of the WCC is driven by the homodimer.

WCC-VVD heterodimerization
As suggested by previous experimental works, the VVD protein is key to the light regulation process through the competitive formation of WCC-VVD heterodimers through interaction of their LOV domains (Chen et al , 2010;Malzahn et al , 2010). In the model, we assume that the light-active forms of the WCC and VVD (WCC* and VVD*) can form a heterodimer (WCC*VVD*) and that the light-active, phosphorylated WCC (WCC* p and VVD*) can also form a heterodimer (WCC* p VVD*).
The lifetime of the heterodimer is governed by: (1) degradation of the complex (l 2 , l 3 ), (2) dissociation of the complex (l −1 ) and (3) degradation of the photoadduct, with the WCC* or VVD* component reverting to their dark form. For (3), the following combinations can occur: VVD dimerizes in a light-dependent manner, forming a rapidly exchanging equilibrium between the monomer and the dimer (Zoltowski and Crane, 2008). Dimerization of VVD competes with heterodimerization, thereby allowing formation of the transcriptionally-active WCC homodimer. We model this in the same way as light-activation of the WCC and homodimerization of the activated WCC, by a two-step process that involves first activation of the VVD protein and then homodimerization of the activated species.

FRQ-mediated phosphorylation of the WCC
The WCC-FRQ negative feedback loop forms the core of the circadian clock. Phosphorylation of the WCC by FRQ inactivates the WCC and therefore, its own transcription.
Experimental evidence has shown that the FRQ level is much lower than that of the WCC (Schafmeier et al , 2005(Schafmeier et al , , 2006 leading to the conclusion that FRQ modulates phosphorylation of WCC, rather than binding to WCC in a 1:1 stoichiometry.

Model equations
Based on the above details, we developed a system of 14 ordinary differential equations describing changes in protein and mRNA concentrations of WCC, VVD and FRQ. The variables are listed in Table S2. Note that we denote WCC by W and VVD by V. An asterix indicates the light-activated form of the protein.

Parameter estimation
To estimate the parameters of the model, we fitted the light-induction data ( Fig. 1) using two different methods: a Bayesian technique and a frequentist inference approach.

Bayesian inference and MCMC
In the Bayesian approach, the parameters are assumed to be random variables that follow a particular distribution. Inferences are made based on the posterior distribution of the parameters, given the data and prior information about the parameters. We first construct the posterior probability distribution where p(·|·) denotes a conditional probability, d is the data and q the parameters to be determined. The left hand of Eq. (1) is the quantity in which we are interested.
The right hand side of Eq.
(1) are quantities that we can calculate: The prior, p(q) summarizes any prior information we have on the possible values of the parameters.
Here we use an uninformative prior for each of the parameters.

p(d|q):
We assume that the error at each data point is Gaussian distributed to give where the d i are the experimental values at each ith time andd i are the estimated values obtained from the set of parameters, q.
A Markov chain with equilibrium distribution p(q|d) is then generated using the Metropolis-Hastings algorithm (Metropolis et al , 1953;Hastings, 1970). New parameter values, q ′ , are randomly drawn from the uniform distribution U (q − δ, q + δ), where δ is the size of the random walk. The new candidate parameters are accepted with probability If the sampling procedure converges for a particular parameter, then the parameter is well-determined by the data and statistics such as the mean and variance of the parameter distribution can be obtained.
We adopted a two-step procedure by first fitting the vvd SS692 light-induction data (Fig. 1). We fixed the fitted parameters and determined the remaining parameters involving VVD by fitting to the wild type data. A subset of parameters was fixed within the range of biochemically reasonable values. Hill coefficients were also restricted. The other parameters for which no such estimates were available were fitted. For our fits, 1 ×10 4 iterations were performed with a burn-in period of usually half the total number of iterations. Samples were collected every tenth iteration thereafter. Parameter values are listed in Tables S3-S7 with their standard deviations. These parameters have been used to simulate the results in the main paper. (1) Max. rate wc-1 transcription 308 hr −1 -K 11 (1) Half-maximal rate wc-1 transcription 1.37 au 0.016 n 11 (1) Hill coefficient wc-1 transcription 0.97 0.0080 k 12 Basal wc-1 transcription 13.49 hr −1 0.19 k 17 wc-1 degradation 8.2 hr −1 - (1) Max. rate frq transcription 250 hr −1 -K 13 (1) Half-maximal rate frq transcription 0.09 a.u n 13 (1) Hill coefficient frq transcription 1.56 k 14 (1),(2) Max.rate frq transcription 56.0 hr −1 -K 14 (2) Half-maximal rate frq transcription 14.61 au 0.21 K 14 ′ (1) Half-maximal rate frq transcription 0.0026 au 0.000016 n 14 (1),(2) Hill coefficient frq transcription 1.92 k 18 frq degradation 9.8 hr −1 (2) W-driven transcription.

Maximum likelihood estimation and profile likelihood analysis
A more traditional approach to parameter fitting is that of maximum likelihood estimation. Here, one first describes an objective function which defines the agreement of the experimental data with the predictions of the model. Commonly, this is the weighted sum of squared residuals: where y D kl denotes d-data points for each variable k at time points t l . The σ D kl are the corresponding measurement errors. The model predictions for a set of parameters θ are given by y k (θ, t l ). The parameters are then estimated byθ giving a point estimate of the parameters that give the best fit of the model to the data. For normally distributed noise, this is also the negative log likelihood of the data being observed for a given set of parameters θ: Given that there are usually a large number of parameters that are to be determined by a much smaller number of experiments, it is important to infer how well the parameters can be identified. The problem of identifiability can be addressed by use of the profile likelihood which can be used to distinguish between structural and practically non-identifiable parameters as well as calculating confidence intervals for the parameters (Raue et al , 2009(Raue et al , , 2011Steiert et al , 2012). The profile likelihood is calculated for each parameter by which means χ 2 (θ) is re-optimized with respect to all remaining parameters θ j =i for each parameter θ i in turn. The shape of the profile likelihood then distinguishes between the different types of identifiability. A threshold, ∆ α can be derived from the log-likelihood of the best fit and the α quantile of the χ 2 distribution to define a confidence interval. A confidence interval of parameter θ i , [σ − i , σ + i ], to a confidence level α implies that the true value of θ i lies within this interval with probability α. A parameter is defined as identifiable if the confidence interval [σ − i , σ + i ] is finite. For a parameter which is identifiable, the χ 2 stays below the threshold ∆ α for a desired confidence level α but then exceeds this threshold at finite values (Raue et al , 2009). If χ 2 (θ) is constant then the parameter is said to exhibit structural non-identifiabiity and its value cannot be determined. Practical non-identifiability occurs if, in one direction, the increase in χ 2 (θ) stays below the threshold ∆ α (Raue et al , 2009). In this case, the confidence interval of the parameter can be finitely bounded in one direction and in the other, it is infinite. An algorithm that can be used as a starting point for calculating the profile likelihood is given in Raue et al (2009). Using this framework, we re-fitted the light-induction data (Fig S4).  Figure S4: Fitting of the light adaptation data using a maximum likelihood optimization approach Fitting shows that the key features of the experimental results can be reproduced (blue curves: wild type; red curves: vvd-mutant). The wild type maintains its responsiveness to the subsequent increases in light intensity, while the successive peaks of the mRNA levels decrease in the vvd-mutant.
The profile likelihoods for the parameters related to VVD (q 6 , q 7 , q 8 ), heterodimerization (l 1 , l −1 ) and the rate of replenishment of the inactive WCC pool, l 5 are shown in Fig. S5. The dashed red line indicates the confidence interval for a confidence level of α = 0.95. The parameters could not be fully identified, although upper and lower bounds could be obtained for some of the parameters. Parameter values of the best fit are listed with corresponding upper or lower limits on their values found from the profile likelihood analysis in the Tables S8-S12. The fitted values using the MCMC algorithm lie within these limits. The profile likelihoods for the parameters related to VVD (q6, q7, q8), heterodimerization (l1, l−1) and the rate of replenishment of the inactive WCC pool, l5. Red circles indicate the best fit. The dashed line indicates the 95% confidence intervals based on the χ 2 distribution; described in Raue et al (2009). Of particular interest is that the rate of replenishment, l 5 , is bounded below and the best fit is found at the upper range of the parameter space explored. The parameter space for l 5 as well as other photoadduct decay variables (l 4 , l 6 and l 7 ) was truncated above at 0.8 hr −1 . This was done because better fits were found with increasing values of l 5 , giving photoadduct decay times of less than an hour, while experimental measurements suggest that the rate of photoadduct decay occurs on the order of hours (Zoltowski et al , 2009 (1) Max. rate wc-1 transcription 49.55 hr −1 >10 K 11 (1) Half-maximal rate wc-1 transcription 0.060 au n 11 (1) Hill coefficient wc-1 transcription 2.4 k 12 Basal wc-1 transcription 10.31 hr −1 <5 k 17 wc-1 degradation 6.42 hr −1 >1 (1) Half-maximal rate frq transcription 0.29 a.u n 13 (1) Hill coefficient frq transcription 1.5 k 14 (1),(2) Max.rate frq transcription 68.35 hr −1 >5 K 14 (2) Half-maximal rate frq transcription 78.20 au - Half-maximal rate frq transcription 0.0010 au n 14 (1),(2) Hill coefficient frq transcription 1.00 k 18 frq degradation 7.44 hr −1 >1 Table S9: Transcription module parameters for FRQ. (1) W*W*-driven transcription.

Supplementary Text 2 Futile cycling provides the means for repeated induction
In order to demonstrate that futile cycling provides a mechanism that could account for repeated sensitivity to increasing step changes in stimuli levels, we constructed a simplified model that represents the core features of adaptation in the Neurospora model (Fig. S6A). Here the inactive species X is activated by an external stimulus to give X*, the activated species which in turn synthesizes its inhibitor I. The activated species is able to form a complex with the the inhibitor. This complex C, (corresponding to the WCC-VVD heterodimer in Neurospora) is able to either dissociate to its two components or dissociate with the activated component X* reverting to its inactivated state, X. The parameter k repl describes the rate of replenishment. We first fitted the vvd mRNA data to this model and showed that this reduced model contains the elements necessary to reproduce the wild type and vvd-mutant behavior: adaptation is achieved in both the wild type and mutant case, with the mutant peaks decreasing in increasing stimulus levels (Fig. S6). (A): Simplified network to highlight the role of futile cycling. Upon stimulation, X is activated to state X*, which then synthesizes the inhibiting species, I. X* and I interact to form the complex C. C can reversibly dissociate but can also dissociate such that the active X* reverts back to its inactive form X. This constitutes the futile cycling component. The parameters ksynx and k degx (marked red) are fixed at 10 a.u hr −1 and 0.01 hr −1 respectively. (B)-(D): The key features of adaptation and the differences between the wild type and the vvd-mutant are reproduced by the reduced model.
We performed a profile likelihood analysis and found that similar to fitting the full Neurospora model, not all parameters can be uniquely determined. However, in contrast, the parameter values are restricted either on an upper or a lower boundary (Fig. S7). Of note is that the parameter crucial to the futile cycle, k repl , is clearly identifiable (Fig. S7A). If the value k repl lies to the right of the confidence interval, then the replenishment rate is too fast and the replenished X is reactivated and therefore, effective adaptation cannot be achieved. On the other hand, if it is too slow, then the wild type tends to the vvd-mutant behavior with diminishing responses to increased stimuli levels. To generalize the result that futile cycling is a method that can be used to achieve repeated responsiveness, we optimized the parameters of the model (Fig. S6) by defining an objective function expressing the qualitative features of adaptation, rather than fitting directly to data. The characteristics of the peak of the responses of a wild type and a mutant system were used to place constraints on the parameters, giving us a constrained optimization problem.
The two key features of adaptation -the ability of a system to respond to step inputs and then to return to the pre stimulus level in the continued presence of the stimulus -can be defined as the sensitivity and precision of a network (Ma et al , 2009). The features that are required to determine the sensitivity and precision are the peak height P i , the adapted level, O i and the input levels, I i (Fig. S8).   Figure S8: Adaptation and repeated sensitivity.
An input I is stepped up successively and at each step increase, the system should respond giving a significant peak, Pi, but then return to its pre-stimulus level, Oi in the continued presence of the input.
Sensitivity is the relative change in the system peak (P 2 ) compared to the steady-state value before the stimulus (O 1 ) with reference to the change in input levels I 1 and I 2 : Precision is given by the difference between the pre-stimulus and post-stimulus steady-state values (defined as the inverse of the relative error by Ma et al (2009)): In contrast to other studies of adaptation, we are interested in the repeated response of a system to further step increases in the stimulus level. Therefore we phrase the sensitivity and precision with respect to n = 4 increasing stimulus levels (Fig. S8). For the objective function, we used the definition of precision: where the precision is defined with respect to the pre-stimulus steady-state value, O 1 .
A key difference between the wild type and the vvd-mutant lies in the magnitude of the responses elicited at each successive, increasing stimulus level. We utilize this fact to define the following constraints to the optimization problem. In the case of the wild type, successive peaks heights should not be diminished: where we set P wildtype = 0.8. In the case of the mutant behavior, the system should no longer be sensitive and the successive peak heights should decrease: where P mutant = 0.1. For both constraints we specify P 3 and P 4 with respect to the first peak, P 2 . We used the Matlab function fmincon with the interior-point algorithm to determine the parameters that solves: subject to the constraints (7) and (8) for the network shown in Fig. S6, where F is the relative error (Eq. (4)) The input levels are: I 1 = 10 −6 ; I 2 = 1; I 3 = 10 and I 4 = 100. The optimization was run for 300 sets of randomly-selected initial parameter values in the interval [10 −4 10]. Similar to constructing a profile likelihood by recalculation of the χ 2 , we also investigated each parameter in turn, starting from the parameter set of the best fit and re-optimizing the system subject to the constraints. While fitting to data provides a natural measure of a good fit (for example, the χ 2 ), here we utilize the fact that the constraints, defining a wild type and vvd-mutant network, must be met and thus sets the criteria for a successful fit. The best fit of the 300 parameter estimation runs is marked by a red circle, while the dashed, vertical lines indicate the boundary at which the constraints are no longer met (Fig. S9). Panels A-C show the parameters related to the futile cycle motif: replenishment of X, and the formation and dissociation of the complex C. While there is no significant change in the minimum value over the range of the parameter search, the rate of replenishment k repl , is bounded below giving a critical value where the constraints are no longer met. The same is also true of the synthesis rate of the inhibitor I, k syni (D). The degradation rate of C, k degc (F), is also a factor in ensuring the constraints are met, as is the degradation of the activated species, X* (G). The constraints specifically distinguish between a wild type and a mutant situation. The boundaries of these constraints show that these parameters must be non-zero, and support the idea that the futile cycling motif is capable of producing adaptation and maintaining responsiveness.