Dynamical modelling of the heat shock response in Chlamydomonas reinhardtii

Global warming is exposing plants to more frequent heat stress, with consequent crop yield reduction. Organisms exposed to large temperature increases protect themselves typically with a heat shock response (HSR). To study the HSR in photosynthetic organisms we present here a data driven mathematical model describing the dynamics of the HSR in the model organism Chlamydomonas reinhartii. Temperature variations are sensed by the accumulation of unfolded proteins, which activates the synthesis of heat shock proteins (HSP) mediated by the heat shock transcription factor HSF1. Our dynamical model employs a system of ordinary differential equations mostly based on mass-action kinetics to study the time evolution of the involved species. The signalling network is inferred from data in the literature, and the multiple experimental data-sets available are used to calibrate the model, which allows to reproduce their qualitative behaviour. With this model we show the ability of the system to adapt to temperatures higher than usual during heat shocks longer than three hours by shifting to a new steady state. We study how the steady state concentrations depend on the temperature at which the steady state is reached. We systematically investigate how the accumulation of HSPs depends on the combination of temperature and duration of the heat shock. We finally investigate the system response to a smooth variation in temperature simulating a hot day.

The green microalgae Chlamydomonas reinhardtii 10 is a widely studied, easy to grow photosynthetic 11 model organism, for which techniques have been 12 whether the interacting molecular mechanisms that 48 have been identified experimentally are sufficient 49 to reproduce and thus explain observed responses.

35
As argued in [Matuszyńska and Ebenhöh, 2015], 36 every mathematical model is usually constructed 37 with a particular purpose and research questions 38 it should help to answer. A model serves as an in 39 silico workbench that can be used to simulate the 40 behaviour of the modelled system in very diverse 41 situations. This allows exploring a variety of sce-42 narios potentially difficult to test, or not yet tested, 43 experimentally. Such simulations often allow the 44 proposition of new hypothesis on the functioning 45 of the biological system under investigation, and to 46 suggest which new experiments could be performed 47 to test these hypotheses. Thus, models provide a 48 theoretical framework to complement the experi-49 mental effort. 50 plast [Willmund and Schroda, 2005]. 1 We develop our mathematical model with the 2 goal to provide a generic description of the heat 3 shock response. Therefore, we decide to simplify 4 the model by including only one generic heat shock 5 protein (indicated by HP), which can represent any 6 HSP present in C. reinhardtii. The simulation re-7 sults are consequently compared with data relative 8 to various of the above mentioned HSPs. In C. rein- 9 hardtii the only HSF (among the two encoded in 10 the genome) known to be activated by heat shock is 11 HSF1, characterized in [Schulz-Raffelt et al., 2007]. 12 Therefore, HSF1 corresponds to the HSF described 13 in our model. In land plants at least 18 different 14 HSF are present [Scharf et al., 2012], which is an 15 example of the fact that in general gene families 16 in C. reinhardtii are smaller than in land plants 17 (see for instance [Schroda, 2004] for another such 18 example concerning chaperones). This simplicity 19 further supports the choice of C. reinhardtii as a 20 suitable model organism to study the heat shock 21 response in plants.  Table 1. In [Schmollinger 27 et al., 2013] it has been shown that the tempera-28 ture increase triggering the HSR is sensed by the 29 accumulation of degenerated proteins P # . Their 30 presence activates a stress kinases SK, which, in 31 the active form SK * , phosphorylates the heat shock 32 factor HSF. The phosphorylated (HSF * ) and un-33 phosphorilated (HSF) heat shock factor can bind 34 to the transcription factor binding sites of various 35 genes, coding for key proteins involved in the HSR, 36 including HSF itself and heat shock protein (HP). 37 In the model, all these genes are described by one 38 variable G, and the transcription of different mR-39 NAs is represented by the individual transcription 40 rates π. Binding of the active form HSF * to genes 41 G induces the production of mRNA coding for heat 42 shock protein (mR HP ) and for the heat shock factor 43 itself (mR F ), whereas the inactive form HSF blocks 44 the transcription. The mRNAs are subsequently 45 translated into the corresponding proteins HP and 46 HSF (with rates η), respectively. The increase in 47 HSF concentration leads to a higher occupation of 48 the gene with the inactive form. The increased con-1 centration of chaperones HP increases the repair of 2 the degenerated protein state P # to their functional 3 form P until a new steady state is reached. 4 The model is described by twelve dynamic vari-5 ables (see Table 1), each representing the concentra-  The equations we employ depend on rate expres-22 sions (see Table 4 of the Supplementary Material) 23 describing the various regulatory processes (ω), acti-24 vation and deactivation steps (ν), synthesis rates (π) 25 and degradation rates (η). For the majority of these 26 rate laws we assume mass action kinetics. We also 27 employ some additional non-linear terms having a 28 behaviour of the type of Michaelis-Menten kinetics 29 or Hill kinetics, listed in Table 4 of the Supplemen-  Heat shock protein (HSP)  Figure 1: Scheme of the signalling network that we use to model the HSR. Temperature T acts via the Arrhenius law ω TP on the protein level P. Higher temperature increases ν P leading to more degenerated proteins P # . This activates stress kinases SK by a hill kinetics ω PS which increases phosphorilization of the heat shock factor HSF. If HSF is bound to the gene G, mRNA for the heat shock factor HSF and for the heat shock protein HP is generated by the corresponding production rates π. The mRNA is translated into the proteins HP and HSF and degraded by rates η.
the heat shock factor itself (panel C) and of the 23 heat shock protein (panel F The simulation results (left panel of Fig. 3) shows 28 the effect of decreasing the rate constant k F on the 29 dynamics of the HSF mRNA concentration. Clearly, 30 decreasing the rate constant leads to a reduced max-31 imal mRNA concentration and a delayed response. 32 As can be seen in the right panel of Fig. 3, the 33 same behaviour is observed in the experiments: an 34 increased inhibitor concentration leads to a reduced 35 and delayed transcription of HSF mRNA. More-36 over, it can be observed that for a small decrease 37 of the rate constant the response is qualitatively 38 unaltered, whereas a larger decrease corresponding 39 to higher Staurosporine concentration has a qualita-40 tive impact, leading to a long delay and considerable 41 reduction of mRNA levels.  Figure 2: Typical behaviour of the model illustrated inducing a HSR via an increase of the temperature from 25°C to 42°C applied at t = 20 min (represented by a red background in the figures). A: Due to temperature increase at t = 20 min functional proteins P are mis-folded leading to an increased P # level. B: The degenerated proteins bring inactive stress kinases SK into their active form SK * . C: Due to active stress kinases, the heat shock factor (HSF) is phosphorylated (HSF * ). D: The heat shock factor HSF binds to free gene loci G, the bound form HSF * G activates mRNA production, and HSF un-binding blocks transcription. E: The initiated gene transcription leads to mRNA production of the HSF and the heat shock protein as shown. F: Due to translation of the corresponding mRNA, the HP concentration increases until the response is switched of. The small degeneration rate of the chaperon leads to a slow decrease after the onset of the HSR. The normalization factors used to represent the concentrations in arbitrary units can be found in

12
Radicicol is a specific inhibitor of HSP90 activity.  which (for wider purposes than studying the HS 10 response) an ARS gene, encoding for the enzyme 11 arylsulfatase (ARS from now on) is placed under 12 control of the HSP70A promoter. The authors show 13 that whenever the HSP70A gene is activated also the 14 ARS enzyme is produced. Under the assumption of 15 a direct proportionality between the concentration 16 of the ARS enzyme and its activity, the authors 17 could monitor the activity of the HSP70A promoter 18 by measuring the ARS activity. This construct was 19 then used to systematically expose C. reinhardtii 20 cells to two subsequent heat shocks to find out what 21 is the minimum time one needs to wait after the 22 first heat shock to observe a full HSR also in the 23 second heat shock. It turned out that a waiting 24 time of around 5 h is required to re-establish the 25 capability to induce a full HSR.    low. However, even this low level is sufficient to 1 induce expression of the HSP genes, so that on the 2 mRNA level the second response appears as equally 3 strong as the first.  To find plausible and realistic parameters, we 10 initially considered biologically reasonable ranges for 11 all rate constants and manually tuned these until we  the experimental data demonstrates that the model 10 is able to reproduce the qualitative, and often the 11 quantitative (relative, not absolute on which we 12 do not have nformation) behaviour of the data ex-13 tremely well (see Sections 3.1 and 3.2).

15
As shown above, our mathematical model, which 16 has been calibrated to control experiments only, 17 can reasonably well reproduce drug treatments as 18 well as the double heat shock experiments. The 19 agreement of simulation results and experimental 20 data therefore supports the notion that our current 21 understanding of the heat shock response is basically 22 correct. The mathematical model can therefore 23 serve as a theoretical framework in which data can 24 be interpreted in a sophisticated and quantitative 25 way. Another purpose of model building is the 26 ability to make novel predictions. We have therefore 27 employed our model to simulate scenarios which give 28 insight into our understanding of the heat shock 29 response, but which are either difficult to test or 30 have not yet been tested experimentally.  32 We first investigate which response the model pre-33 dicts upon exposure to a prolonged HS, and how 34 the system adapts to persistently high tempera-35 tures. Experimentally, the systems-wide response 36 to long-term HS was investigated in [Hemme et al., 37 2014], where cells adapted to 25°C were exposed to 38 42°C for a period of 24 h, followed by 8 h at 25°C 39 (recovery phase).

40
The simulation results for this scenario, where 41 the temperature increase was simulated at time 42 t = 0, are shown in Fig. 6. Two distinct phases of 43 the response can be distinguished. The first, early 44 HS phase, lasting approximately 3 h after applying 45 the heat shock, represents the initial heat shock 46   not too high, remarkably the concentration of un-10 folded proteins at steady state is kept very low by 11 the response at the level of all the other species, very 12 close to zero and in particular well below one per-13 cent of the total amount of proteins. On the other 14 hand, when the temperature increases considerably 15 the HSR is no more able to efficiently counteract 16 the accumulation of degenerated proteins which 17 accumulates at concentrations high enough to kill 18 the cell. This accumulation is evident in panel A 19 of Fig. 7, a magnification of which is provided in 20 Fig. 8. Finally, we have also verified that, for each 21 of the values of the temperature that we have con-22 sidered, the model exhibits a realistic stationary 23 behaviour, i.e. the associated steady state (which 24 is a non-equilibrium one) is stable. To do so, firstly 25 we computed numerically the Jacobian of the vector 26 field associated to the ODEs system summarized in 27 Table 3 of the Supplementary Material. Next, we 28 evaluated the Jacobian of the system at the steady 29 state under investigation. Then we computed the 30 eigenvalues of the Jacobian to determine the stabil-31 ity of the steady state. We repeated the procedure 32 for the steady state associated to each of the values 33 of the temperature considered. We obtained that all 34 the nine eigenvalues have always negative real part, 35 showing that the steady state is always a stable one. 36

38
The questions which we want to address in this 39 section are the following: "is C. reinhardtii more 40 stressed for a short HS with high temperature or 41 a long one but with a lower temperature? What 42 is the trade-off between the two?". A related ques-43 tion is "does the production of HSP occurs only 44 under very intense HS conditions (high tempera-45 ture, long duration) or does it occurs also for very 46 small temperature increases or very short HS?".   To reproduce the typical experimental set-up em-34 ployed in many studies we have considered through-35 out this work heat shocks provided by means of a 36 step-wise increase of the temperature from a lower 37 value to an upper value. These are conditions that 38 can be easily reproduced experimentally, and allows 39 to make straightforward tests in vivo. They could 40 also be employed for possible applications of HS (as 41 for instance expressing proteins of interest. Never-42 theless, these are situations not so likely to occur 43 often in nature.

44
Thus, which kind of HS is a C. reinhardtii cell 45 going to experience in the wild? This green algae 46 is widely distributed around the world in various 47 environments such as soil and fresh water. Thus, 48

19
In Section 4 we have used the model to simulate 20 situations not yet tested experimentally, from which 21 we derive interesting results. We have shown in 22 Section 4.1 that the system can adapt to higher 23 temperatures during heat shocks longer than three 24 hours, by shifting to a new steady state. Two dis-25 tinct phases are clearly visible: an early HS lasting 26 for about the first 3 h, and a late HS in which the 27 system shows adaptation (a new steady state is 28 reached). The recovery phase is characterized by a 29 recovery of the conditions pre-HS that occurs over 30 several hours. 31 We have studied in Section 4.2 the variation of 32 the steady state concentrations w.r.t. changes in 33 the temperature. The concentrations of HP, as 34 well as of the mRNAs, increase with increasing 35 temperature, but for not too high temperatures the 36 concentration of unfolded proteins P # is kept very 37 close to zero, in particular well below 1% of the total 38 amount of proteins [P ] + P # . For temperatures 39 too high the HSR cannot prevent the accumulation 40 of degenerated proteins and the cell dies. 41 We have used the model to systematically investi-42 gate how the accumulation of HSPs depends on the 43 combination of temperature and duration of the heat 44 shock, in Section 4.3. Short (smaller than 10 min) 45 HSs do not provide enough time for a significant 46 response at the level of [HP ], the maximum of [HP ] 47 for any given temperature is obtained at around 80 48 to 100 min, after that a somewhat smaller [HP ] is 49 reached and maintained, and for long enough HS a 50 higher temperature provides higher [HP ]. 1 We have finally investigated the system response          Let us remark that the values of the variables are initiated at the initial conditions above, but before applying any HS we let the system run for a long time, so that it has reached the steady state when we apply any HS. This part of each simulation is not shown in the plots.

Rate law
Reaction Table 4: Kinetic rate laws used in the model. The reactions are those represented in the scheme of Fig. 1, and the rate laws are mainly based on mass action kinetics, a part from some terms winch follows Arrhenius law or have a Michaelis-Menten or Hill kinetics behaviour. The last is ω P S = [P # ] m P m 0 +[P # ] m , where m = 5 and P 0 = 600 µM. The term following the Arrhenius law is ω T P = A · exp − Ea R·T , with an activation energy of E a = 174.440 KJ mol −1 , perfectly in the wide range reported in the literature for the activation energies of protein denaturation due to thermal stress (as discussed in [Bischof andHe, 2006,He andBischof, 2003]), and with A = 9.4318 × 10 28 , R = 8.3144598 J mol −1 K −1 the ideal gas constant and T the temperature. Finally, the basal rates are computed as p RF basal = d RF · d F · 0.02125 µM/k π F and p RP basal = d RP · d HP · 17.5 µM/k π HP .  In this section we present some additional simula- is necessary to allow the system to quickly refold 10 unfolded proteins. The HSR corresponding to a 11 second HS occurring at 3 h 30 min after the first 12 is similar, but enhanced. When the second HS 13 occurs 5 h or more after the first we see that the 14 concentration of HSP is approaching the level it had 15 before the HS (panel F), all the other quantities are 16 approximately back to the original values, and an 17 almost full HSR takes now place when the second 18 HS is applied (panels B, C, D and E).

19
It is very interesting to remark that, while the 20 concentrations of all the species go back to the values 21 that they had before the first HS quite fast after 22 the end of the first HS, the HSP does not (panel 23 F), and this allows to avoid during the second HS 24 having any but a tiny amount of unfolded protein 25 P # with respect to the amount during the first HS. 26 This can be seen in panel A of any of the sub-figures 27 of Fig. 12, where the concentration of degenerated 28 proteins P # increases by a considerable amount 29 during the first HS, while considerably less during 30 the second even when this is occurring many hours 31 after the first.

32
The behaviour we observe in our simulations likely 33 indicates that the production of HSF which follows 34 a first HSR and the accumulation and slow degrada-35 tion of HSP have the role of preparing the system 36 for a subsequent occurrence of the same stressing 37 situation (HS) already encountered in the past, thus 38 representing a transient molecular memory.   Fig. 6b). This com-11 parison requires the simplifying, but still reasonable, 12 additional assumption that the activity of the ARS 13 enzyme is roughly directly proportional to its con-14 centration. Given the simplicity of the extension 15 that we used to include this into our model, and the 16 rough estimation of the involved parameters, the 17 qualitative agreement between simulation and data 18 is already remarkable.

C Supplementary Material:
ex-20 tended description of the calibra-21 tion of the model 22 In this section we provide more details on the proce-23 dure that we have employed to calibrate the model. 24 We describe in detail the objective function used and 25 which data are involved in calibration. We then de-26 scribe the random sampling of the parameter space 27 and the gradient search employed to determine the 28 final set of parameter values.  Table 5).   Very low values of the RMS can occur everywhere 10 in the interval used for the random scan, depending 11 on the values assumed by the other parameters. 12 We subsequently plotted the values of the RMS 13 as function of each couple of parameters (obtaining 14 400 figures). We found that sometimes there is 15 some correlation or anti-correlations between the 16 preferred values for the two parameters of the couple 17 (an example is shown in Fig. 15). Nevertheless, for 18 the majority of the couples no such correlation can 19 be identified.  Fig. 5). Defining such an RMS is much more 11 arbitrary than defining the RMS w.r.t. the controls 12 of the feeding experiments, because of the previ-13 ously mentioned hypothesis on the proportionality 14 between the enzyme concentration and its activity. 15 We do not combine the two RMSs, as this would 16 require to attribute to the two a weight which would 17 be highly arbitrary.

18
For these reasons we only compute this second 19 RMS a posteriori as a check, for the best 5000 20 among the random parameters sets selected above, 21 and show the distribution of the values of the two 22 different RMS in Fig .17. This shows that minimiz-23 ing both at the same time cannot be obtained, and 24 one needs to find a trade-off between the two. The 25  for the 20 rate constants employed in our model (see 6   Table 2). We compute the corresponding value of 7 the root mean square RM S ( x 0 ). We then compute 8 numerically the gradient of the RMS at that point 9 ∇RM S ( x 0 ) (by approximating partial derivatives 10 using the symmetric difference quotient). 11 We then proceed along the direction op-12 posite to the gradient towards a new point 13 x n+1 = x n − γ ∇RM S ( x n ) in the parameter space 14 which provides a smaller value of the RMS. We do 15 so iteratively until a termination criterion described 16 above is satisfied, and label the iteration number by 17 n. At each iteration, we need to decide which is the 18 length of the step γ that we want to use in the direc-19 tion opposite to the gradient. To do so, we imple-20 ment a line search with the aim to loosely minimize 21 the function f (γ) .
w.r.t. γ, i.e. along the direction opposite to the 23 gradient. This means finding the value of γ which 24 minimizes the function f (γ). We do so numerically 25 employing a modification of the bisection rule based 26 on the Golden ratio to save computation time.

27
Since the orders of magnitude of the parameters 28 are very different, we expect the isosurfaces of the 29 function RM S ( x) to be far away from being spher-30 ical. This would lead to a very slow convergence of 31 double HS, as can be seen from Fig. 18 and Fig. 19.   Fig. 21 and discussed in the Supplementary 12 Material.

13
The figures mentioned above show the correspond-14 ing results of the simulations performed with our 15 model using the final parameter set determined as 16 described in the previous section. The experimental 17 set-ups, the salient features of each experiment and 18 the way to use our model to simulate these exper-19 iments have been widely described in Sections 3.1 20 and 3.2. As we can see from the figures, the model 21 reproduces well the qualitative and the (relative) 22 quantitative behaviour of these experimental data. 23

D Supplementary Material: Maxi-24
mal accumulation of degenerated 25 proteins as function of time neces-26 sary to increase the temperature 27 during HS 28 We might wander how the maximum concentration 29 of unfolded proteins accumulated during a HS de-30 pends on how fast the temperature has increased 31 from the initial value T low to the final value T high . 32 To study this we simulate what happens when pro-33 viding to the system a HS by acclimating the system 34 at the temperature T low =25°C, then starting from 35 t = t 0 , then increasing the temperature following a 36 cosinusoidal function until when the temperature 37 T high =42°C is reached at t = t 0 + τ , and then 38 keeping the temperature at T high . In this way the 39 function T (t) is not only continuous, but also every-40 where differentiable. We repeat for various values 41 of the time τ which is the time required for the 42 temperature to increase from T low to T high , and we 43 show the result in Fig. 20. We observe that for in-44 stantaneous increase in temperature, up to increases 45 which requires about 1 minute, the maximum value 46  the absolute values, we nevertheless know given its 10 definition that the z-score is linear with respect to 11 the concentration value represented by each data 12 point, then we do not expect any distortion on the 13 vertical axes of the plot if we would be able to 14 transfer these data to the corresponding original 15 values, which justify the comparison of Fig. 21. The 16 qualitative agreement is good, nevertheless we can 17 notice a somewhat faster increase during the first 18 100 minutes of the model's prediction w.r.t. the 19 data.