Using a model comparison approach to describe the assembly pathway for histone H1

Histones H1 or linker histones are highly dynamic proteins that diffuse throughout the cell nucleus and associate with chromatin (DNA and associated proteins). This binding interaction of histone H1 with the chromatin is thought to regulate chromatin organization and DNA accessibility to transcription factors and has been proven to involve a kinetic process characterized by a population that associates weakly with chromatin and rapidly dissociates and another population that resides at a binding site for up to several minutes before dissociating. When considering differences between these two classes of interactions in a mathematical model for the purpose of describing and quantifying the dynamics of histone H1, it becomes apparent that there could be several assembly pathways that explain the kinetic data obtained in living cells. In this work, we model these different pathways using systems of reaction-diffusion equations and carry out a model comparison analysis using FRAP (fluorescence recovery after photobleaching) experimental data from different histone H1 variants to determine the most feasible mechanism to explain histone H1 binding to chromatin. The analysis favors four different chromatin assembly pathways for histone H1 which share common features and provide meaningful biological information on histone H1 dynamics. We show, using perturbation analysis, that the explicit consideration of high- and low-affinity associations of histone H1 with chromatin in the favored assembly pathways improves the interpretation of histone H1 experimental FRAP data. To illustrate the results, we use one of the favored models to assess the kinetic changes of histone H1 after core histone hyperacetylation, and conclude that this post-transcriptional modification does not affect significantly the transition of histone H1 from a weakly bound state to a tightly bound state.

S3 File: Perturbation analysis for models 6, 7, and 8 Two-population model The system of reaction-diffusion equations describing the two population model (model 0) is given by where D eff is the effective diffusion coefficient, and k on and k off are the binding and unbinding rates, respectively.

Model 6
The system of reaction-diffusion equations describing model 6 is given by (2) Assuming that γ b = λ b ε , and γ u = λu ε , where ε << 1, the reaction-diffusion system (2) can be rewritten as If we now consider a perturbation expansion for u(x, t), w(x, t) and v(x, t) of the form we obtain the following leading-order system for (3) From the second equation in (5), we obtain the quasi-steady state relation w 0 ( , and use the quasi-steady state, we note that w 0 (x, t) = γ 1+γ c 0 and u 0 (x, t) = 1 1+γ c 0 . Substituting this into the first equation, the leading-order system (5) becomes the reaction-diffusion system of two equations This is equivalent to (1) with Thus, we conclude that if the turnover of weakly bound biomolecules to/from a freely diffusing state is sufficiently fast the reaction-diffusion system of three equations (2) can be approximated with the reactiondiffusion system of two equations (6). If we use the solution of the reaction-diffusion system of equations (2) to fit histone H1.5 data assuming a diffusion coefficient D = 25 µm 2 /s, we obtain parameter values of If we use those values to estimate parameters in (1), we obtain D eff = 3.227462, κ on = 0.001011, κ off = 0.002231.
In Fig. 1 (left panel) we show the recovery curve according to both model 6 (2) with parameter values (7), and its leading-order approximation (6) equivalent to the two population model (1) with parameter values (8).
Alternativetly, we could fit first the leading order approximation (6), and then estimate parameters in model 6 (5) according to the relation However, parameters γ b and γ u cannot be uniquely determined. By fitting directly the leading order approximation (6) to the same histone H1.5 data, we obtain the following parameter estimates Assuming a diffusion coefficient D = 25 µm 2 /s, we can estimate the following parameters for model 6 γ = 5242.045, η b = 0.00021, η u = 0.00112.
In Fig. 1 (right panel), we used these parameter estimates to sketch the explicit solution of the full model 6 for different values of γ u and γ b , keeping their ratio γ fixed. three-population model 6 (blue) and its two population leading-order approximation (green) to histone H1.5 data. Right panel: Recovery curves of the three population model 6 approximating H1.5 data for several parameter estimates with values greater than two orders of magnitude.

Model 7
The system of reaction-diffusion equations describing model 7 is given by Assuming that γ b = λ b ε , and γ u = λu ε , where ε << 1, the reaction-diffusion system (11) can be rewritten as If we now consider a perturbation expansion for u(x, t), w(x, t) and v(x, t) of the form we obtain the following leading-order system for (12) From the second equation in (14), we obtain the quasi-steady state relation w 0 (x, t) = γu 0 (x, t), where γ = λ b λu = γ b γu . If we define c 0 (x, t) = u 0 (x, t) + w 0 (x, t), and use the quasi-steady state, we note that w 0 (x, t) = γ 1+γ c 0 and u 0 (x, t) = 1 1+γ c 0 . Substituting this into the first equation, the leading-order system (14) becomes the reaction-diffusion system of two equations This is equivalent to (1) with Thus, we conclude that if the turnover of weakly bound biomolecules to/from a freely diffusing state is sufficiently fast the reaction-diffusion system of three equations (11) can be approximated with the reactiondiffusion system of two equations (15). If we use the solution of the reaction-diffusion system of equations (11) to fit histone H1.5 data assuming a diffusion coefficient D = 25 µm 2 /s, we obtain parameter values of κ b = 0.009542, κ u = 0.002045, γ b = 0.090338, γ u = 0.016153 If we use those values to estimate parameters in (1), we obtain D eff = 3.792223, κ on = 0.001447, κ off = 0.002045.
In Fig. 2 (left panel) we show the recovery curve according to both model 7 (11) with parameter values (16), and its leading-order approximation (15) equivalent to the two population model (1) with parameter values (17). Alternativetly, we could fit first the leading order approximation (15), and then estimate parameters in model 7 (14) according to the relation However, parameters γ b and γ u cannot be uniquely determined. By fitting directly the leading order approximation (15) to the same histone H1.5 data, we obtain the following parameter estimates D eff = 0.0047682216, κ on = 0.0002122792, κ off = 0.0011231743.

Model 8
The system of reaction-diffusion equations describing model 8 is given by Assuming that γ b = λ b ε , and γ u = λu ε , where ε << 1, the reaction-diffusion system (20) can be rewritten as If we now consider a perturbation expansion for u(x, t), w(x, t) and v(x, t) of the form we obtain the following leading-order system for (21) From the second equation in (23), we obtain the quasi-steady state relation w 0 ( , and use the quasi-steady state, we note that w 0 (x, t) = γ 1+γ c 0 and u 0 (x, t) = 1 1+γ c 0 . Substituting this into the first equation, the leading-order system (23) becomes the reaction-diffusion system of two equations This is equivalent to (1) with Thus, we conclude that if the turnover of weakly bound biomolecules to/from a freely diffusing state is sufficiently fast the reaction-diffusion system of three equations (20) can be approximated with the reactiondiffusion system of two equations (24). If we use the solution of the reaction-diffusion system of equations (20) to fit histone H1.5 data assuming a diffusion coefficient D = 25 µm 2 /s, we obtain parameter values of If we use those values to estimate parameters in (1), we obtain D eff = 3.749684, κ on = 0.001249, κ off = 0.002045.
In Fig. 3 However, parameters γ b and γ u cannot be uniquely determined. By fitting directly the leading order approximation (24) to the same histone H1.5 data, we obtain the following parameter estimates D eff = 0.0047682216, κ on = 0.0002122792, κ off = 0.0011231743.
In Fig. 3 (right panel), we used these parameter estimates to sketch the explicit solution of the full model 8 for different values of γ u and γ b , keeping their ratio γ fixed.

Remarks
As explained in the manuscript, Figs. 1-3 show that the recovery curve obtained from fitting any of the threepopulation models is significantly different from their leading-order approximations for parameter estimates not greater than two orders of magnitude, but equivalent for parameter estimates greater two orders of magnitude. This occurs because the turnover of weakly bound biomolecules to/from a freely diffusing state is not sufficiently fast when the rapid interaction is characterized with parameters values not greater than two orders of magnitude; in other words, the assumption of high turnover rates (γ b and γ u ) is not met.