BiP Binding to the ER-Stress Sensor Ire1 Tunes the Homeostatic Behavior of the Unfolded Protein Response

Computational modeling and experimentation in the unfolded protein response reveals a role for the ER-resident chaperone protein BiP in fine-tuning the system's response dynamics.


Model Variables
The variables in the model describe the number of molecules rather than concentrations. For ease of notation and clarity of model equations, we will put square brackets around each variable, e.g. [X]. A description of the model variables is given in the

Model Assumptions
In building the computational model, we adopted the following biologically motivated assumptions -Ire1 is activated through contact with unfolded protein that are not in the process of folding properly, that is U , U d , and U d · B. It is not activated by U · B which is folding properly.
-DTT breaks disulfide bonds in the unfolded protein and folding complexes and this process occurs as an enzymatic reaction. We also assume that there is a constant pool of DTT diffusing in and out of the ER, replenishing the local population.
-Ero1 mends broken disulfide bonds resulting from DTT. The population of unfolded proteins with broken disulfide bonds vs. those without is dependent on the relative populations of Ero1 vs. DTT. In reality, Ero1 is one protein of many that assist in the formation of disulfide bonds, others include oxidoreductases from the PDI family. For simplicity, Ero1 will represent the entire mechanism for disulfide bond formation and will act as an enzyme.
-Ero1 and DTT are equally efficient in their enzymatic reactions. Therefore the competition is 1-to-1.
-Unfolded proteins regulate an allosteric switch in Ire1 from the inactive state to the active state. Ire1 then allosterically deactivates based on the active state's deactivation rate. This model is chosen over the alternate model that unfolded protein activate Ire1 by binding to it and deactivate by dissociating from it. Mathematically, the allosteric model is simpler and because of the small number of Ire1 and the generally large number of unfolded proteins, either scenario would yield essentially the same results.
-The deactivation rate of active Ire1 is non-linearly dependent on the population of active Ire1. The relationship we use is pictorially described in Figure S10a. Ire1 molecules are known to form foci which are correlated with Hac1 mRNA splicing. The nonlinear function we employ in the model is meant to serve as a simple phenomenological description of the deactivation rate of Ire1 from the Ire1 foci, formed during stress. This relationship can describe a variety of underlying mechanistic schemes of Ire1 inactivation -The population of Ire1 is approximately 256 [Ghaemmaghami et al. (2003)].
-We are neglecting the lectin chaperone pathway.
-We are neglecting the ATP/ADP cycle of BiP which is known to switch BiP's affinity (high/ATP and low/ADP) to molecules such as Ire1 and unfolded proteins. In this model we set BiP's affinity constant and equal to both Ire1 and unfolded proteins.
-We are neglecting the contribution of ERAD.

Chemical Reactions and Reaction Rate Constants Used in Model
The reaction rate constants we will use in the equations are the traditional reaction rate constant divided by the volume of the space involved (area if on a membrane), e.g. our c µ = k µ /V where k µ is the rate constant for some species µ. The model consists of three modules. Below, we provide a description of the chemical reactions modeled and the reactions rates used.

Protein Folding Module
Reaction Description BiP attachment to unfolded protein to form folding complex.
Completion of folding and release of folded protein from folding complex.
Breakage of unfolded protein's disulfide bonds by DTT.
Formation of disulfide bonds by Ero1 in the unfolded protein with broken disulfide bonds.
BiP attachment to unfolded protein with broken disulfide bonds to form misfolding complex.
Breakage of disulfide bonds in folding complex by DTT.  [Berg and Purcell(1977)]. D c = 1 µm 2 s −1 is a typical cytosolic diffusion coefficient. a p = .028 µm is the typical protein size in the ER. V er = 2.15 µm 3 is the approximate volume of the ER.  [Axelsen and Sneppen(2004)].

Ire1 activation, Hac1 splicing, and splicing reporter module
Reaction Description BiP binding to Ire1 to form inactive complex.
Splicing of Hac1 mRNA and reporter mRNA [I0] n [I 0 ] n +[I 1A ] n represents a non-linear (active Ire1 cooperative) decay rate as illustrated in Figure S10b. We have set γ [I 1A ] ≈ γ [I 1 ·B] and n = 4.5, I 0 = 45. For more details see Crucial Parameter Fitting section below.
β [H u γ [H 1 m ] -8.33 × 10 −4 s −1 . Decay rate of Hac1 mRNA corresponding to a mean decay time of 20 minutes as determined from [Sidrauski et al.(1996)] (Figure 5c and page 408). We use the same value for both spliced and unspliced.
Splicing rate of Hac1 mRNA corresponding to an 11 minute splicing time inferred from Figure 3 of [Mori et al.(1997)], which also shows that the total number of mHac1, spliced and unspliced remains relatively constant, an assumption in our model. The total number of molecules being spliced per second is is required since the minimum of these quantities will determine the number of molecules being instantaneously spliced.
Decay rate of the splicing reporter mRNA. This is assumed to be the same as γ [H 1 m ] .

BiP and Ero1 transcription module
Reaction Description is the transcription hill function for the UPRE promoter where a 0 = 296.5 and a 1 = 5.26. The function is fitted (see Figure S8) to the normalized data in [Credle et al.(2005)  γ [B] -1.39 × 10 −4 s −1 . Decay rate of BiP corresponding to a 2 hour mean decay time referenced from the modeling paper [Axelsen and Sneppen(2004)].

Kinetic Model Equations
Module 1: protein folding dynamics Module 2: Ire1,Hac1 mRNA, and splicing reporter dynamics Module 3: BiP and Ero1 dynamics To simulate the various mutants, the equations were adjusted as follows: Starting from the equilibrium solution for non-stressed conditions (zero DTT), the DTT levels (represented by D 0 as the molecular number of DTT within the ER) were adjusted as a function of time, simultaneously solving the equations using the ordinary differential equation solver, ODE15s, in MATLAB.

Crucial Fitted Parameters
In the section, we discuss parameters for which we did not find any quantitative measurements or inferences in the literature. There are seven such parameters. To fit these parameters, we use a subset of the data describing the dose dependent splicing reporter experimental data sampled at 200 minutes after stress induction for the wild type, hac1∆, ire1 bipless , and hac1∆ + ire1 bipless . Our best fit is shown in Figure S9a. Figure S9b shows the results for the modeled washout experiment based on our fitted model predicting a delay in the ire1 bipless mutant, which was also experimentally verified. The details for parameter fitting are listed below.
S u -The source rate of unfolding proteins, whose value is important relative to the value of total BiP. The best fitted results were obtained when the basal level of total BiP (≈ 430, 000) was set to be 15-20 percent above the basal level of folding proteins. This value results in the population of free BiP to be around 60, 000. In the model, the larger the amount of free BiP to bound BiP, the slower the onset of splicing in the wild type system in the 5 mM DTT experiment. This is due to the fact free BiP sequesters Ire1 until there is enough unfolded proteins in the ER to titrate free BiP. Therefore, the choice of this parameter was set to match well with the onset of max splicing (20-30 minutes). However, in principle, this value can be changed without affecting the qualitative behavior of the system, especially the delay observed in the washout experiment. Instead, different values for S u affect the magnitude of the delay. Figure S11a illustrates the sensitivity of the model to this parameter by varying the value of S u . Here we plot the modeled washout experiment for S u = 370 (dashed), 310 (solid), and 250 (dashed). All values produce a delayed response in the ire1 bipless case relative to the wild-type, but the delay dependent on the value of S u . [B] ), was fit to data from hac1∆ and hac1∆+ ire1 bipless . This parameter was especially picked to fit the half maximal induction of the dose response curves for .4 and .25 mM DTT. Since these experiments have no transcriptional feedback, the basal level reaches approximately 1,350,000 molecules per ER volume, keeping in mind that Ero1 in this model represents the full system of molecules combating the effects of a similar number of DTT molecules.
N [B m ] -N [B m ] = 4 was found to best fit the dose response data for the transcription feedback response of the wild type and ire1 bipless experiments. This value is within the range of the measured transcriptional reporters.
N [E 1 m ] -N [E 1 m ] = 7 was found to best fit the dose response data for the transcription feedback response of the wild type and ire1 bipless . This value is within the range that transcriptional reporters measured.
The constants c [I·B] , γ [I·B] , c up and γ [I·U ] dictate the system's dynamics, and determine the levels of free BiP and misfolded/unfolded proteins that can compete for free Ire1.
for the unfolded protein and c [I·B] [B] for BiP. We found that the ratio α up = c [I·B] /c up = 150 works best in fitting the model to data. This value makes intuitive sense since both Ire1 and the unfolded protein (even those with BiP attached) will generally diffuse much slower than free BiP. Figure S11b shows that the quantitative aspects of the model are dependent on the variation to α up for the values 100 (dashed), 150 (solid), and 250 (dotted). Here the intuition is that the lower the value of α up , the lower the size of the relative population of misfolded proteins that is needed to compete for free Ire1.
γ [I 1 ·B] -γ [I 1 ·B] represents the dissociation rate of BiP from the inactive complex. This quantitiy is important since it determines the amount of time that an Ire1 molecule is free versus bound to BiP in a given condition of ER stress. Although no real quantitative measurements have been recorded, there are some estimates based on the work in [Bertolotti et al.(2000)] where BiP association to Ire1α in mammalian cells was studied. We will assume similar result for yeast.
In the experiments in [Bertolotti et al.(2000)] for non-stressed cells, Ire1α/BiP complex and Ire1α were immunoprecipitated. It was found that the majority of Ire1α were in the complex form with BiP. Based on this, we can assume that at equilibrium, the ratio of the "ON" to "OFF" rates are the same as the ratio of the number of Ire1/BiP to free Ire1, neglecting active = R.
In our model we set R = 12.5. Figure S11c plots the results of the washout experiments for 3 values of R: 125 (dashed), 12.5 (solid), 1.25 (dotted). In this sensitivity test, for a given value of R, we also scaled γ [U ·B] and γ [I 1A ] by the relative change in R from its best fit value of 12.5, i.e 10 (a), 1 (b), .1 (c), respectively. The results in Figure S11c show that the system is robust for values of R ≥ 12.5, but at R = 1.25, the full splicing response to 5 mM DTT does not occur. An explanation for this is that the basal level of active Ire1 (splicing) are higher in the non-stressed case for R = 1.25, thus are preadapted for higher levels of stress.  Figure S10b. The most important assumption of Ire1 deactivation that was crucial for fitting the model to the data was that the active state of Ire1 was stabilized as a function of increasing the active Ire1 population. For our model we were able to fit the data by setting γ [I 1A ] = γ [I 1 ·B] , with a best fit of n = 4.5, I 0 = 45. Reasonable fits of the model were obtained over the ranges of n = 3 − 6, I 0 = 20 − 60. However, for n = 1, I 0 = 15, Figure S9c-d plots the results of the dose dependent splicing reporter sampled at 200 minutes as well the results for the washout experiment. Both the switch-like behavior the ire1 bipless dose response and the delay in the ire1 bipless washout experiment are lost. That said, we would like to point out that this cooperative function is a simplified representation of the ire1 cooperativity in the foci dynamics.

Additional Parameter Sensitivity Analysis
The modeling analysis implemented in the paper assumed the mean protein folding time, 1/γ f old , to be 20 minutes. However, the ratio S u /γ f old dictates the basal unfolded protein population. Therefore, we checked that the qualitative behavior of the model is insensitive to the folding time by varying the mean between 12 and 20 minutes, while keeping the ratio (basal population) constant. Figures S12a-b show the dose response curves at 200 minutes and the washout experiments, respectively, for folding times of 12, 15 and 20 (best fit) minutes. In all cases, the system exhibited the same qualitative trends seen in the data. We also examined the sensitivity of the model to different values of the cytosolic diffusion coefficient, D c , which affects many reaction constants simultaneously. We varied D c over a range of .5 < D c < 2.0 µm 2 s −1 . Over this range, the model robustly reproduced the trends seen in the data as seen in Figure S12c for the washout experiments. The dose response curves at 200 minutes (not shown) were all very similar as well.
Our model assumes an Ire1 molecular count of 256 molecules/cell [Ghaemmaghami et al. (2003)], which is based on an average value over many cells. We varied the range of Ire1 from 232 to 281. A notable feature of the model is that while the deactivation dynamics of the wild type system are still faster than the ire1 bipless mutant in the simulation of the washout experiment, a large variability in the shutoff delay was was observed as the number of Ire1 molecules varied for the ire1 bipless mutant ( Figure S12d). This observation constitutes an interesting experimental prediction.