Unearthing the transition rates between photoreceptor conformers

Background Obtaining accurate estimates of biological or enzymatic reaction rates is critical in understanding the design principles of a network and how biological processes can be experimentally manipulated on demand. In many cases experimental limitations mean that some enzymatic rates cannot be measured directly, requiring mathematical algorithms to estimate them. Here, we describe a methodology that calculates rates at which light-regulated proteins switch between conformational states. We focus our analysis on the phytochrome family of photoreceptors found in cyanobacteria, plants and many optogenetic tools. Phytochrome proteins change between active (P A) and inactive (P I) states at rates that are proportional to photoconversion cross-sections and influenced by light quality, light intensity, thermal reactions and dimerisation. This work presents a method that can accurately calculate these photoconversion cross-sections in the presence of multiple non-light regulated reactions. Results Our approach to calculating the photoconversion cross-sections comprises three steps: i) calculate the thermal reversion reaction rate(s); ii) develop search spaces from which all possible sets of photoconversion cross-sections exist, and; iii) estimate extinction coefficients that describe our absorption spectra. We confirm that the presented approach yields accurate results through the use of simulated test cases. Our test cases were further expanded to more realistic scenarios where noise, multiple thermal reactions and dimerisation are considered. Finally, we present the photoconversion cross-sections of an Arabidopsis phyB N-terminal fragment commonly used in optogenetic tools. Conclusions The calculation of photoconversion cross-sections has implications for both photoreceptor and synthetic biologists. Our method allows, for the first time, direct comparisons of photoconversion cross-sections and response speeds of photoreceptors in different cellular environments and synthetic tools. Due to the generality of our procedure, as shown by the application to multiple test cases, the photoconversion cross-sections and quantum yields of any photoreceptor might now, in principle, be obtained. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0368-y) contains supplementary material, which is available to authorized users.


A mathematical model for phytochrome undergoing two thermal reversion processes
In the main text, we have shown a number of examples highlighting the utility of our methodology. Notably, cases whereby a population of photoreceptors contains two sub-populations undergoing different thermal reversion rates. One could envisage that this is due to multiple reactions within a monomeric photoreceptor or that different conformers of photoreceptor dimers undergo different thermal reversion reactions as is the case for full-length phyB from Arabidopsis [1]. Since the case of photoreceptor dimerisation has been treated mathematically elsewhere (see [1] and Section 1.2), we present a model whereby a monomeric photoreceptor undergoes multiple thermal reversion reactions.
The model readsṖ where k 1 and k 2 are the light induced reactions that are proportional to the photoconversion cross-sections (k i = ∑ λ σ λ i z λ λ w as described in [2][3][4] and z λ w is the photon distribution from the light source as described in the main text), β 1 and β 2 are the two thermal reversion rates, γ is the rate of complex association and δ is the rate of dissociation. Note that synthesis and degradation rates have been ignored, since these processes do not occur within our photoreceptor samples that are purified and measured in solution. Here P is assumed to be 'free' photoreceptor molecules, whereas C could take one of two interpretations. First, C is some form of complexed photoreceptor with some unknown interaction partner, H. For example, in our system phyB-N is thought to exist as a monomer. Therefore, H must either be very small or is intrinsically part of the light reaction mechanisms of phyB-N, which could be the case for light-sensing chromophores [5,6]. In our model, C could represent a phyB-chromophore complex that alters the dynamics between conformational states. On the other hand, P and C could be distinct sub-populations of photoreceptor proteins whereby the protein folding machinery produces distinct isomers. For example, dimerization within the protein structure of phyB is known to occur through the GAF domain [7,8]. Thus, P could be one form of phytochrome whereby the GAF domains are unbound, but C represents a sub-population where GAF domains are bound, altering conformational switches. The aim of our work in this section is to solve this system to provide a general form of the equation that describes P A decay in the presence of light and thermal reactions.
The total amount of phytochrome in the system is P 0 = P I + P A + C I + C A and the total amount of H, H 0 = H + C I + C A , is conserved. Therefore, by settinġ X =Ṗ I +Ṗ A andẎ =Ċ I +Ċ A one can find that where 0 ≤ Y ss ≤ P 0 . Using these results we can then re-write our system in terms of P A and Z = P A + C Ȧ We solve the inhomogeneous equatioṅ for u ij = v ij e ψ j t . The resulting eigenvalues and eigenvectors are Using the initial condition the general solution is The solution Z(t) is dependent on ρ ± that, in turn, depends on a number of system parameters, including thermal reversion rates and interaction rates, that can be difficult to elucidate experimentally. Thus, if we assume that two forms of P A naturally exist (as has been postulated experimentally [7,8]), we can set δ = γ = 0 in equation (2) such that ρ + = −(k 1 + k 2 + β 1 ) and ρ − = −(k 1 + k 2 + β 2 ). This is important, since we can estimate the values of β k from our thermal reversion optimisation algorithm allowing us to use our approach to analyse these systems.
All that is left, then, is to calculate the value of the constants A, B, C, and D. Importantly, we wish to know how the percentage of active photoreceptor changes with time. Thus, we can rescale (3) usingx = x/P 0 for Z and the exponential constants. Therefore, we getẐ explained in the main text) and, similarly, is the steady state fraction of P A subpopulation k. From our thermal reversion optimisation routine, as well as obtaining the multiple thermal reversion rates, we also obtain the percentage of the active photoreceptor population undergoing that particular thermal reversion reaction, α 1 and α 2 = 1 − α 1 .
Thus, we can satisfy (4) byÂ Finally, by setting α = α 1 and k i = ∑ λ σ λ i z λ λ w for a given light treatment λ w , we obtain Thus, if α = 1, then B = D = 0 and our system will only describe a single photoreceptor form with thermal reversion rate β 1 . In Figures 3, 4 & 6 of the main text we show the results of our approach given different numbers of thermal reaction equations following (5) or (7).

Forms of steady-state percentages of active photoreceptor, c λ w
In this section we shall present some of the known forms for the steady state percentage of active photoreceptor in a population given illumination wavelength λ w .
As noted in [2], in the absence of thermal reversion, synthesis and degradation, phytochrome dynamics are described by At steady-state, one finds that where k i = ∑ λ σ λ i z λ λ w as before. When the photoreceptor undergoes a single thermal reversion reaction the system reads where β is the thermal reversion rate.
In our study we discuss two further cases that will influence the distribution of active photoreceptors within a population, namely: the existence of two sub-populations of monomeric photoreceptor, and; the dimerisation of photoreceptors. Using the model above in equation (3), and setting t = 0 or t → ∞, yields for the case of multiple photoreceptor sub-populations where β k represents multiple thermal reversion rates. The form of c k λ w can be calculated from our model (equation (1)) at steady state when δ = γ = 0.
To obtain the steady state percentage of active photoreceptor for the case of photoreceptor dimers, we use the model of [1] in the absence of synthesis and degradation reactions. The system reads where P I I is the P I − P I homodimer, P AI is the P A − P I heterodimer, and P AA is the P A − P A homodimer. From these equations we obtain where c 1 λ w is the steady-state percentage of P A -P I heterodimers and c 2 λ w is the steadystate percentage of P A -P A homodimers. Here, the multiplication by n c = 2 of c 2 λ w is required as the homodimer complex contains two light-absorbing chromophores within it's structure and ∑ 2 k=1 c k (σ, Φ, α, β, Ω P , Ω E ) = n c from equation (1) of the main text.
Depending on the case under investigation or what is known about the photoreceptor from previous studies will determine which form of c λ w needs to be used within the algorithm. A further example can be found in Section 1.5.1.

Relating c's to X's in Verméglio method
In the main text, we discuss the construction of search spaces for a photoreceptor's photoconversion cross-sections due to unknown parameters, such as the percentage of (in)active photoreceptor, quantum yields and thermal reversion rates. However, if some of these are known (for example there is no thermal reversion and absorption spectra show that 100% of the photoreceptor is in a particular confirmation) then the values of X 1 and X 2 can be calculated directly. Here, in the case of a photoreceptor with two light regulated species, we show how the X's are related to the steady-state fraction of active photoreceptor after saturating light conditons (c λ w ).
From the definition of absorption spectra and how the photoconversion cross-sections are related to multiple spectra, we have, for any number of available chromophores n c where n cˆ By assuming that experimental spectra share the same relationship to their extinction coefficients as simulated spectra (i.e. that our model is correct) and substituting (9) into (10), one finds that which leads to and so Furthermore, these values satisfy the constraints

Application to a simpler photoreceptor system
As well as the phytochromes, further photoreceptors are studied in plant and optogenetic systems. For example, blue light-and UV-B-responsive photoreceptors have been used in a range of optogenetic tools (see the review in [9]). Importantly, compared to phytochromes, these systems only respond to one region of the light spectrum and can only be de-activated by periods of darkness and thermal relaxation. Consequently, the mathematics can be simplified. Here, we will provide the equations and how to find the optimal solution.
First, since these photoreceptors have only one light-induced state, the equation for absorption can be written as where σ λ I and Φ I are the photoconversion cross-sections and quantum yield of the inactive photoreceptor to reach the active confirmation. Note that here c I = n c − c A with n c = 1 as the photoreceptor exists as a monomer.
Since, after prolonged darkness, the percentage of active photoreceptor is zero and there is no light reaction-induced absorption, we can only obtain the cross section after illumination of activating light conditions. Thus, c λ I = 0 and in the assumed case where a single population of active photoreceptor exist. Furthermore, terms describing the fraction of active photoreceptor after saturating light conditions (t → ∞) are simplified as well. For example, Finally, we note that the Verméglio method is simplified also, sincê (16) which upon substitution into (14) with As with the case of phytochromes, a search space forˆ λ I can be obtained by varying Φ I and calculating all possible values of X 2 . Using this search space, the optimal Φ I and set of σ λ I can be obtained by minimizing the differences between simulated data and measured absorption spectra.

Generalising method to photocycles containing N species
Whilst in the previous section we have simplified the analysis of the main paper for a photoreceptor with a single state, we can also generalise our method to N-step photocycles. Such complex photoreactions have already been shown to govern the dynamics of phytochromes [10,11]. Here, we shall first present the conditions that need to be satisfied to determine the search spaces forˆ λ i and values for R max i (the upper bound for the quantum yield ratios) for i = 1, ..., N before illustrating these results for a 3-step cycle that incorporates the phytochrome species Lumi-R [10][11][12]. Note that here we assume there is no dimerisation between photoreceptors such that n c = ∑ N k=1 c k (σ, Φ, α, β, Ω P , Ω E ) = 1.
From equation (1) of the main text, we see that absorption spectra can be described as a linear combination of extinction coefficients such that whereˆ λ k = 2.303lc tot λ k . Furthermore, by again assuming that our model of absorption is correct such that measured spectra share the same relationship with extinction coefficients as simulated spectra, we can rearrange equation (4) of the main text to read Thus, from equations (17) and (18) we obtain an expression for any steady-state absorption spectra after illumination with light condition z λ where (κA) k is the k-th row of κA. Solving this system for z λ λ A and z λ λ I , and using the relation that c A ∞ + c I ∞ = 1, provides us with the functions of equation (6) in the main text to obtain values of X 1 , X 2 and R max .
Given this reformulation of the problem, one can now generalize the system to a reaction cycle containing N species. In this instance we write where z λ j are the N light conditions required to obtain the N steady state absorption spectra for each species and γ j i are scaling factors such that ∑ j κ kj = 1.
From equation (20) we thus have, for any steady-state absorption spectra, From this, we now obtain sets of simultaneous equations that satisfy this relationship and can be solved to provide us with values of X i and R max i = Φ i /Φ 1 , the ratio of quantum yields for each species relative to the inactive state. For example, In the case of two species, and using c 1 ∞ + c 2 ∞ = 1 we obtain equivalent relations to those of equation (6) in the main text.

Example: the three-species cycle of phytochromes
As an illustration of the N species extension outlined above we can now obtain the relationships for a three species phytochrome system. In this system, inactive P I is photoconverted into the Lumi-R (L I ) state under red light that goes on to form P A . As with previously, P A can then revert back to P I using far-red light or thermally in darkness [10][11][12]. The ODE model of this system is then where k i are light regulated reactions and β j are thermal reactions.
Consequently, one can obtain analytical expression for the steady states of the three componentsP Then, given that (8) of the main text, where σ λ i are the photoconversion cross-sections, Φ i are the quantum yields and λ i are the extinction coefficients of each species i ∈ {I, L, A}, we can rewrite these expressions aŝ Using the conditions in equation (22) we then have which, upon substitution of equation (25) gives Then substituting the quantum yield ratios R A = Φ A /Φ I and R L = Φ L /Φ I and solving the equations one obtains These simultaneous equations can now be solved numerically within the framework of our optimisation procedure to find the optimal quantum yields and photoconversion cross-sections for multiple phytochrome species. However, one should note that to analytically obtain these expressions becomes increasingly challenging as the number of species in the reaction cycle increases [13]. Thus, one may wish to consider other methods of numerically obtaining the percentages of subspecies (possibly through numerical integration of the ODEs given specific parameter values) followed by obtaining numerical solutions of the simultaneous constraint equations (equation (22)).

Calculation of phyB-N photoconversion cross-sections using Butler's method
In order to show the differences between our method and those previously used to estimate the photoconversion cross-sections of full-length phytochrome, we have calculated the quantum yields and photoconversion cross-sections using Butler's method. Essentially, using this method we are assuming that c λ I = 0 (there is no P A after saturating exposure to λ I ), which means that X 1 = 1 from equation (11). The full derivation of the method can be found in [14] but we shall summarise it here.
The first step is to plot L against ln 100E λ w t , where E λ w t is described in equation (2) of the main text given illumination of light with wavelength λ w during the experiment. Here, where z λ is our photon distribution, l 0 is the cuvette width, l is the corrected pathlength, ∆t is the time of illumination, V s is the volume of the sample andĀ λ is the intensity normalised absorption of a blank sample (see Supplementary Table 2 for these parameter values) [14,15].
Upon calculating these functions, the aim is to determine for which value of L allows 100E λ w t = 36.8 (referred to as L 36.8 , for an example see the result for our data in Supplementary Figure 6C) [14,16]. The term L 36.8 refers to the amount of exposure required for absorption to decay to e −1 (∼ 0.368) of the initial value under the assumption that illumination at λ w leads to a single photoproduct (i.e. P A ) [16]. By taking 1 L 36.8 , we get the values σ 657 I = 43.9 m 2 /mol and σ 735 A = 10.5 m 2 /mol using our data calculated at the peaks of photon absorption (656 nm LEDs peak at 657 nm, 740 nm LEDs peak at 735 nm, see Supplementary Figure 2). These are in broad agreement with the scale of values for the photoconversion cross-sections at these wavelengths that we obtained from our algorithm (Supplementary Table 4), espe-cially when compared with the 'Mancinelli' spectra where both our algorithm and Butler's method predict much smaller photoconversion cross-sections for phyB-N (see Figure 6 of the main text). The next step is to perform the following derivations to obtain the quantum yields and a scaling factor between the absorption spectra and photoconversion crosssections.
, where A e (λ A |Ω P , {(t → ∞, z λ λ I )}) is the absorbance of our sample measured at wavelength λ A after saturating inactivating light z λ λ I (and the same for activating conditions).
We obtain a value of X λ A eq = 0.7457 that leads to X 2 ∼ −0.3410.
-below we shall relate this to (16) and show that they are equivalent.
• Φ A = σ λ I A /α λ A . By calculating α λ A one is able to obtain the quantum yield Φ A and then from step 1 calculate Φ I .
Following this procedure, we calculated that Φ I = 0.006 and Φ A = 0.008 (Supplementary Figure 6B), or roughly one fifth of the the values calculated by our optimisation algorithm ( Figure 6 of the main text).
All that remains is to calculate the wavelength-dependent photoconversion crosssections of σ λ A using step 2. This is related to the Verméglio method since ). The resulting spectra can be seen in Supplementary Figure 6A. What is obvious from this figure is that the wavelength-dependency of σ λ A using Butler's method is similar to that obtained from our algorithm (compare the solid and dashed lines of Supplementary Figure 6A), however the relative amplitudes and absolute values at the peak of the photoconversion cross-sections obtained from Butler's method are smaller than those obtained by our optimisation algorithm. Using the calculated set of photoconversion cross-sections and quantum yields, we simulated difference spectra and compared the output to our data and the simulated spectra obtained from our optimal values (Supplementary Figure 6D). From this figure it is clear to see that Butler's model has not calculated values that provide a close match to the time-dependent spectra obtained experimentally and that our optimal values predict.

Absorption spectra optimisation algorithm
Here, we provide pseudocode for our optimisation algorithm that aims to minimise the difference between simulated and measured absorption spectra for given sets of photoconversion cross-sections.
1. Define the upper and lower bound of the estimates. We have used: 2, N] in the N species case) 2. Start optimization routine to minimize the difference between simulated and experimental absorption spectra. We used the 'Interior-point' algorithm of the f mincon optimization function in MATLAB.
3. Define the initial conditions. To do this we used the rng and rand functions within MATLAB to generate random values between the upper and lower bounds for each parameter.
a. set Φ A = RΦ I (or Φ i = R i Φ 1 in the N species case).
b. calculate the values of X 1 and X 2 (or X i ) for given quantum yields (see note below).
c. interpolate search space of photoconversion cross-sections for these values of X's. d. calculate the cost score as defined by equation (9) of the main text.
e. If Ω is less than the previous estimate, keep it as the current best score and try the next parameter estimates.
f. Once Ω reaches pre-defined termination criteria, stop the optimization and return the minimal score. In our work we have set the termination criteria to be: 6. If the algorithm is boot-strapped further, then re-start points 4-7.
From the resulting optimal parameter sets the mean σ λ -spectra can be obtained by calculating X 1 and X 2 for the optimal Φ's and interpolating the search space of photoconversion cross-sections.
Note: As stated in the main text, the functions of X 1 and X 2 are implicit. This means that they need to solved numerically. In this work we have used the MATLAB numerical solver vpasolve. However, in some unforeseen cases, the numerical solver fails to converge to an accurate solution of X 1 and X 2 for a given set of Φ. At these points, we obtain numerical solutions in the local vicinity of the original Φ values (perturbing the values by 0.0001 in each direction of the Φ-space) and then interpolate these values to obtain a solution of X 1 and X 2 at the original values of Φ. This, in effect, produces a continuous space of potential X 1 and X 2 values. Due to the small size of the perturbations, we believe that this process does not greatly impact our results.

Thermal reversion optimisation algorithm
The algorithm to calculate thermal reversion parameters works as follows: 1. Use equation (3) of the main text to set A D t = ∑ n k=1 α k e −β k t D . Define the number of n exponential functions that should be fit.
2. Define the upper and lower bound of our parameter estimates. We have used: β k ∈ 0.0001, 0.01 . For cases where a single thermal reversion rate is required, α can be fixed to 1 or β values can be fixed to zero.
3. Perform the optimisation routine 1000 times: a. define initial conditions for the parameter search as random points between the upper and lower bounds; b. use the 'Interior-point' option for f mincon in MATLAB to minimise where E D t is described by equation (2) of the main text and ζ t is the standard deviation of the data (similar to ν λ t above). c. once Ω reaches pre-defined termination criteria, stop the optimization and return the minimal score. In our work we have set the termination criteria to be: i. maximum number of iterations ( MaxIter ) = 10000.

Protein analysis protocol
Final preparations of phyB-N containing the N-terminal 650 amino acids were stored in 50 mM sodium phosphate buffer containing 300 mM NaCl, pH 7.0. The concentration of phyB was determined using BCA assay and purity was analysed according to SDS-PAGE. Absorption measurements were performed using Agilent/Hewlett-Packard 8453 Diode-Array UV-Vis Spectrophotometer (Agilent Technologies, CA, USA).

Absorption spectra measurements of photoswitching
Photoswitching experiments were conducted by diluting phyB-N in a buffer (50 mM sodium phosphate buffer containing 300 mM NaCl, pH 7.0) to final concentration 0.34 mg/ml as determined from BCA analysis. This concentration is below an optical density of 0.2 in order to avoid inner filter effects and ensuring homogeneous illumination of the sample [20]. Absorption and illumination measurements were performed in 500 µl Hellma Quartz (Hellma GmbH & Co. KG, Müllheim, Germany) black absorption cuvettes that are 1cm thick. Illumination was performed either using 656 nm or 740 nm light emitting diodes (LEDs) at 20 cm distance yielding 11 µmol/m 2 s flux rates, respectively. The photon distribution of the LEDs was determined using a Jobin Yvon-Spex Fluorolog 3.22 (Horiba Scientific, NJ, USA) in combination with micromax fibre optics. The head of the fibre was placed 1.5 m in front of the LED and spectra were recorded. Recording of red to far-red experiments were performed by placing the cuvet with phyB-N under 656 nm LEDs for 5 min followed by 740 nm illumination for different time periods. Absorption spectra were recorded after 10, 20, 30, 40, 50, 60, 80, 100, 120, 160, 200 s of 740 nm illumination. The far-red to red experiment was conducted similarly, starting with 8 min of 740 nm light followed by 5,10,15,20,25,30,40,50,60, 120 and 180 s periods of 656 nm illumination. Measurements were conducted at room temperature (∼22 • C) and recorded in biological duplicate.

Thermal reversion experiments
After 5 min illumination of 656 nm light at 22 • C, a cuvet containing phyB-N was placed in complete darkness within the spectrophotometer. Absorption spectra were obtained on a binary logarithmic scale starting at 42 seconds until a final measurement after 172,032 s (∼48 h). Measurements were conducted at 22 • C in duplicate.

Supplementary Tables
Supplementary Table 1  V is the estimated sample volume for each simulated experiment chosen such that output of Butler's model matches input parameters t λ x = time period for illumination of wavelength λ X