The Zimm-Bragg model

1The main assignment¶

The Zimm-Bragg model is a simple model to describe the cooperative coil-to-helix transition that occurs in polypeptides. The theoretical background, discussed during an earlier lecture, can be found here, and a Jupyter notebook containing a code that solves the model (theoretically, as well as numerically through Monte Carlo simulations) is here. Here you are asked to adapt my code (or write your own) to study the effect of some of the model parameters on the melting temperature of the transition, defined as the temperature at which half of the residues are in the helical state, i.e. $\theta = 0.5$.

In order to be positively evaluated, your code should be able to produce plots (or the data, to be plotted separately) of the melting temperature as a function of σ, $N$, $\Delta H$, and $\Delta S$ (where the latter are the entalphic and entropic contributions to $s$). The accompanying paper (or README) should also contain an interpretation of the results. My advice is to use the exact solution rather than Monte Carlo simulations.

Note that the melting temperature should be plotted against the quantity we wish to vary, while keeping all the others fixed to some sensible value.

2Possible extensions¶

• Compare the results with MC simulations. Within the statistical noise, do you obtain the same results? If not, comment why. Hint: the larger $N$, the more steps you need. The same is true if the temperature is low, or if the other parameters are such that the Boltzmann factor tends to be very small, since in these cases the probability to accept a flip becomes very small, slowing the sampling of indepedent configurations, and therefore worsening the quality of ensemble averages.
• As discussed during the lectures, the propensity to form secondary structures, and helices in particular, varies from amino acid to amino acid. This property can be incorporated in the model by assigning different values of $s$ (and therefore $\Delta H$ and $\Delta S$) to different residues, depending on their identity. An interesting exercise in this regard would be:
  • Extend the code to support residues of two (or more, if you are brave) types, using $s_1$ and $s_2 \neq s_1$ as their elongation parameters.
  • Run MC simulations and compare the estimated value of θ at fixed temperature and some sensible values of the other parameters, with that obtained with the original code, simulated with $s = \frac{s_1 + s_2}{2}$.
  • Here we have introduced disorder in the system: you will have to average over different sequences to compute the value of θ that will be used for the comparison. This procedure is called averaging over the quenched disorder, and it is a widespread concept in advanced statistical mechanics.

3Additional details¶

The melting temperature can be estimated by evaluating the melting curve (as done in the notebook), and then use it to find the temperature at which the curve crosses the $\theta = 0.5$ line. Specifically, since our $\theta = \theta(T)$ curve is monotonic^[1], then we can invert the function to obtain $T = T(\theta)$, and then compute the melting temperature as $T_M = T(0.5)$. However, in order to do so we have to “transform” the discrete data points into a continuous function. This operation is called “interpolation”, and there are plenty of libraries available to do this in any programming language. Here I will show how the interpolation can be carried out with scipy, which is a widely-used Python library.

We first import the interp1d function:

from scipy.interpolate import interp1d

Calling it with the $x$ and $y$ coordinates of our data points used as parameters will return the interpolated function, $T(\theta)$:

interpolated_theta = interp1d(theta, Ts)

Now we can call the function to evaluate at what temperature there is the given fraction of helical residues:

print("The melting temperature is", interpolated_theta(0.5))

Here is a full working example:

from scipy.interpolate import interp1d

# test data: they must have the same length
Ts = (0.2, 0.4, 0.6, 0.8, 1.0, 1.2)
theta = (1.0, 0.92, 0.61, 0.1, 0.01, 0)

interpolated_theta = interp1d(Ts, theta)
print("The melting temperature is", interpolated_theta(0.5))