Minimum energy paths for conformational changes of viral capsids

TL;DR

This paper applies Large Deviations Theory to model the conformational changes of a viral capsid, identifying the most probable transition paths and analyzing the stability of intermediate states in a simplified dodecahedral model.

Contribution

It introduces a novel application of Large Deviations Theory to study energy paths in viral capsid conformational changes, focusing on a simplified geometric model.

Findings

Computed minimum energy paths for capsid transition

Identified stable intermediate states during conformational change

Provided insights into the stability of weakly interacting protein structures

Abstract

In this work we study how a viral capsid can change conformation using techniques of Large Deviations Theory for stochastic differential equations. The viral capsid is a model of a complex system in which many units - the proteins forming the capsomers - interact by weak forces to form a structure with exceptional mechanical resistance. The destabilization of such a structure is interesting both per se, since it is related either to infection or maturation processes, and because it yields insights into the stability of complex structures in which the constitutive elements interact by weak attractive forces. We focus here on a simplified model of a dodecahederal viral capsid, and assume that the capsomers are rigid plaquettes with one degree of freedom each. We compute the most probable transition path from the closed capsid to the final configuration using minimum energy paths, and discuss the stability of intermediate states.

Tables

Table 1. Table 1: Half-barriers between states a 𝑎 a to e 9 subscript 𝑒 9 e_{9} .

	$b$	$c_1$	$c_2$	$c_3$	$d_1$	$d_2$	$d_3$	$d_4$	$d_5$	$e_1$	$e_2$	$e_3$	$e_4$	$e_5$	$e_6$	$e_7$	$e_8$	$e_9$
$a$	0.61
$b$		0.42	0.60	0.61
$c_1$					0.25	0.42	0.60	0.59
$c_2$						0.24	0.42	0.41	0.59
$c_3$							0.41
$d_1$										0.25	0.41				0.59
$d_2$										0.09	0.24	0.41	0.59	0.57
$d_3$												0.23	0.40		0.24			0.41
$d_4$											0.08			0.40	0.25	0.57	0.41
$d_5$													0.23	0.22			0.40

Table 2. Table 2: Half-barriers between states f 1 subscript 𝑓 1 f_{1} to j 𝑗 j .

	$f_1$	$f_2$	$f_3$	$f_4$	$f_5$	$f_6$	$f_7$	$g_1$	$g_2$	$g_3$	$g_4$	$g_5$	$h$	$i$	$j$
$e_1$	0.24		0.57			0.41
$e_2$	0.09			0.56		0.24
$e_3$			0.39			0.08
$e_4$		0.22	0.07	0.22			0.36
$e_5$			0.09				0.37			0.0484
$e_6$			0.23	0.39	0.38	0.08
$e_7$				0.07			0.20
$e_8$					0.22		0.37		0.0475
$e_9$		0.21			0.20
$f_1$								0.09	0.24	0.23		0.55
$f_2$											0.07		0.02
$f_3$												0.21	0.05
$f_4$											0.07	0.06
$f_5$														0.01
$f_6$									0.08		0.38
$f_7$															0.002
$g_1$													0.23
$g_2$													0.08
$g_3$													0.08
$g_4$															0.02
$g_5$															0.02
$h$														0.07
$i$															0.05

Equations

$$A_{ij}=\left\{\begin{array}[]{ll}1&\text{if }ij\in E\\ 0&\text{otherwise},\end{array}\right.\qquad i,j\in\{1,\dots,12\}.$$

$$\rho (g){\mathbf{e}}_j:={\mathbf{e}}_{\sigma (g)(j)},g\in \mathcal{I},$$

$$\tilde{\rho }=\rho \otimes \Gamma ,$$

$$\hat{\rho }=\underset{i=1,2,3,5}{⨁}{\rho }_i\otimes \Gamma ,$$

$$E(\mathbf{x})=\sum _{i=1}^{12}f(x_i)+\gamma \sum _{i,j=1}^{12}{A}_{ij}g(x_i,x_j),\gamma \in \mathbb{R},$$

$$E(H\mathbf{x})=E(\mathbf{x}),\forall H\in {\mathcal{H}}_3,\forall \mathbf{x}\in {\mathbb{R}}^{12},$$

$$\nabla (E(H\mathbf{x}))=H\nabla E(\mathbf{x}),\forall H\in {\mathcal{H}}_3.$$

$$\frac{\partial E}{\partial x_i}(\mathbf{x})=f^{\prime }(x_i)+2\gamma \sum _{j=1}^{12}{A}_{ij}{g}_1(x_i,x_j),i=1,\dots ,12.$$

$$\begin{array}{rl}& \frac{{\partial }^{2}E}{\partial x_i^2}(\mathbf{x})=f^{\prime \prime }(x_i)+2\gamma \sum _{j=1}^{12}{A}_{ij}{g}_{11}(x_i,x_j),\\ & \frac{{\partial }^{2}E}{\partial x_i\partial x_j}(\mathbf{x})=2\gamma A_{ij}{g}_{12}(x_i,x_j).\end{array}$$

$$h(x)=f^{\prime }(x)+10\gamma g_1(x,x)=0.$$

$$aI+bA,$$

$$h^{\prime }(x)=a(x)+5b(x)={\mu }_1,$$

$$\dot{\mathbf{x}}=-\nabla E(\mathbf{x}),$$

$$d\mathbf{x}=-\nabla E(\mathbf{x})dt+ϵd\mathbf{W},$$

$$S_T(\mathbf{\phi })={\int }_0^T|\dot{\mathbf{\phi }}+\nabla E(\mathbf{\phi }){|}^{2}dt,$$

$$\underset{ϵ\to 0}{\lim }{ϵ}^2\ln {\mathbb{P}}^{\mathbf{x}}(\mathbf{x}(T)\in D)=-\underset{\mathbf{\phi }\in {\mathcal{C}}_{\mathbf{x},D}}{\min }{S}_T(\mathbf{\phi }),$$

$${\mathbf{\phi }}^{\ast }=\underset{\mathbf{\phi }\in {\mathcal{C}}_{\mathbf{x},D}}{\mathrm{argmin}}{S}_T(\mathbf{\phi }),$$

$$\lim_{\epsilon\to 0}{\mathbb{P}}^{{\boldsymbol{x}}}\left(\sup_{0\leq t\leq T}|{\boldsymbol{\varphi}}^{*}(t)-{\boldsymbol{x}}(t)|<\delta\Big{|}{\boldsymbol{x}}(T)\in D\right)=1.$$

$$V(\mathbf{x},\mathbf{y})=\underset{T>0}{\inf }\underset{\mathbf{\phi }\in {\mathcal{C}}_{\mathbf{x},\mathbf{y}}}{\min }{S}_T(\mathbf{\phi }),$$

$$p(\mathbf{x})\asymp \frac{1}{Z}\exp \left(-\frac{V({\mathbf{x}}_0,\mathbf{x})}{{ϵ}^2}\right),$$

$$V({\mathbf{x}}_0,\mathbf{x})=2(E(\mathbf{x})-E({\mathbf{x}}_0)).$$

$$(\nabla E{)}^{\perp }({\mathbf{\phi }}^{\ast })=0,$$

$$(\nabla E{)}^{\perp }({\mathbf{\phi }}^{\ast })=\nabla E({\mathbf{\phi }}^{\ast })-⟨\nabla E({\mathbf{\phi }}^{\ast }),\hat{\mathbf{\tau }}⟩\hat{\mathbf{\tau }},$$

$$Q_{\alpha\beta}=\left\{\begin{array}[]{ll}\displaystyle\exp\left(-\frac{2}{\epsilon^{2}}(E({\boldsymbol{x}}_{\alpha\beta})-E({\boldsymbol{x}}_{\alpha}))\right)&\text{if ${\boldsymbol{x}}_{\alpha}$ and ${\boldsymbol{x}}_{\beta}$ are connected}\\ \displaystyle 0&\text{otherwise}\end{array}\right.$$

$$Q_{\alpha \alpha }=-\sum _{\beta \ne \alpha }{Q}_{\alpha \beta }.$$

$${\Pi }_{\alpha \alpha }=0,{\Pi }_{\alpha \beta }=-\frac{Q_{\alpha \beta }}{Q_{\alpha \alpha }},\alpha \ne \beta .$$

$${\Pi }_0=\underset{ϵ\to 0}{\lim }\Pi .$$

$${\boldsymbol{s}}\in S\Leftrightarrow\left\{\begin{array}[]{l}\displaystyle\lim_{t\to+\infty}{\boldsymbol{x}}(t,{\boldsymbol{s}}))=\hat{\boldsymbol{x}}({\boldsymbol{s}}),\\ s_{i}=1\quad\Leftrightarrow\quad\hat{x}_{i}({\boldsymbol{s}})=1,\\ s_{i}=0\quad\Leftrightarrow\quad\hat{x}_{i}({\boldsymbol{s}})<1,\end{array}\right.$$

$$\beta (\mathbf{s},{\mathbf{s}}^{\prime })=V(\hat{\mathbf{x}}(\mathbf{s}),{\mathbf{x}}_m)=2(E({\mathbf{x}}_m)-E(\hat{\mathbf{x}}(\mathbf{s}))),\text{if}\mathbf{s}\to {\mathbf{s}}^{\prime },$$

$$Q({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})=\left\{\begin{array}[]{ll}\exp\left(-\frac{1}{\epsilon^{2}}\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})\right)&\text{if ${\boldsymbol{s}}\to{\boldsymbol{s}}^{\prime}$},\\ 0&\text{otherwise},\end{array}\right.$$

Full text

Structural transitions, viral capsids, large deviation theory

I Introduction

Viral capsids are interesting biological structures assembled from repeated copies of the same protein crick . They are very efficient at their purpose of protecting the genetic material of the virus from the environment, since they are quite stable for a wide range of environmental conditions (cf. zandi_2 for a review).

In order to release the genome inside the host cell, however, the capsid must be able to change configuration and /or disassemble in response to changes in the chemical environment or the interaction with receptors of the host. Such conformational changes often involve the opening of pores in the viral shell through which the nucleic acid exits the virus and is released into the cell (see, for instance, tuthill ; bakker ; blaas2000 ; blaas2002 ; blaas2012 ). More often, though, the capsid is believed to simply disassemble as a consequence of the weakening of the bonds between the capsomers cao .

Also, many viral capsids undergo structural changes during maturation hk97_2 . The assembly of the capsid is often a multi-stage process that may involve various steps towards the infective, final form of the virion. Once the capsomers have assembled to form a closed shell, called procapsid, these still have to undergo conformational changes involving protein cleavage, subunit rotation and /or deformation, and substantial bond disruption and reforming, in order to reach the final, stable form of the infective virus conway .

Hence, the definition of conformational change is necessarily somewhat loose, ranging from the simple mutual detachment of the proteins leading to the complete disgregation of the capsid, to the cleavage of the capsomers triggering complex relative rigid motions of the protein subunits, to the formation of new bonds with changes in the quaternary structure of the proteins.

In any event, the question arises as to which is the basic physics underlying the stability and structural plasticity of the capsid. The forces driving the conformational changes are diverse, and not always known. For instance, the nucleic acid is highly compressed within the closed capsid, and this generates an internal pressure that could destabilize the shell conway , even though in RNA viruses its large negative charge might contribute to the stability of the capsid through Coulomb attraction with the coat proteins. In other cases, there are charges on the faces of adjacent proteins that are masked by ions in the stable, closed capsid but, upon pH changes, lose the ions and trigger a Coulomb repulsion between the proteins speir . Sometimes, instead, it is the cleavage of a particular protein that acts as switch for a conformational change of the capsomers that reach a new stable configuration castellanos .

Cohesive forces that keep the capsid together are also diverse: there are no covalent bonds between the capsomers, and binding usually involves strong hydrophobic bonds between chains ’knotted’ around or inside beta barrels of adjacent proteins, hydrogen bonds, Coulomb attraction, and so on roos .

What is clear, though, is the cooperative nature of the stability of the capsid. The bonds keeping the capsomer together must be comparatively weak in order to allow for easy and fast configurational changes, and the secret of the capsid stability must lie in the joint stabilizing action of adjacent capsomers upon each other bakker ; hk97_2 ; castellanos ; hk97_1 ; Zlotnick2003 ; Zlotnick2009 .

In this work we continue the investigation started in CIT ; giuliana1 (cf. also (bruinsma2015 )), to study under which conditions configurational changes in viral capsids involve either simultaneous collective movements, or completely disordered unrelated events, or a cascade of local destabilization events leading to a wavefront propagating along the shell. Indeed, the energy cascade hypothesis has been proposed in ross and experimentally demonstrated in cardone for the maturation of HK97, a complex process that involves multiple (possibly icosahedral) intermediates.

We employ a very simple model of a dodecahedral viral capsid, in which the capsomers are pentagonal rigid units endowed with a single degree of freedom, that may represent a geometric variable, such as a displacement or a rotation of the capsomer as in ERAV and CCMV, or an internal variable measuring the conformation of the protein subunits as in HRV or HK97. Due to its simplicity, the dodecahedral capsid model is widely used in the literature to perform simulations, but its use here is motivated by the ERAV capsid bakker , a pseudo $T=3$ capsid in which disassembly does occur by the relative motion of almost rigid pentagonal units made of 20 coat proteins. Our approach could also be applied to larger, more complex capsids, in which both the structure and the protein-protein interactions are richer. The purpose here is to put the energy cascade hypothesis on a firm theoretical ground in the simplest setting as possible.

We assume that the configurational change is the result of the competition between a driving force, that that we view either as the internal pressure due to the confinement of the genetic material inside the capsid, or the Coulomb repulsion between charges at the capsomer-capsomer interfaces, and a counteracting interaction term describing the attraction between adjacent capsomers, consistent with the form of capsid energy proposed in kegel2004 for HBV. Notwithstanding the extreme simplicity of the interaction forces, the geometry of the interactions between adjacent capsomers is responsible for the high complexity of the energy landscape, that has many minima corresponding to metastable states of the capsid.

As mentioned above, conformational changes can be triggered by a number of diverse factors. In this work we focus on transitions driven by changes in the chemical enviroment of the capsid that modulate the relative intensity of the two main energetic contributions to capsid stability: for instance, the removal of ions due to changes of the pH usually unmasks charges that are responsible for augmented Coulomb repulsion between the interfaces. All these factors are reckoned by a single parameter $\gamma$ in the energy that weighs the relative strength of the competing forces, and is such that when it is above a certain threshold, the original configuration of the capsid is stable (immature provirion in maturation events or closed capsid in the disgregation case), and when it decreases below this threshold the initial configuration is unstable. We show that when the critical threshold is reached, the transition between the initial and the final state occurs directly, without intermediates, and requires a large energy CIT . Hence, we assume that the parameter $\gamma$ is near but above the critical threshold, in order to retain the fine features of the inter-capsomer interactions in our model.

In order to explore the energy landscape and determine the minimum energy paths between metastable states, we use here the formalism of Large Deviations Theory freidlin . The formalism provides both a solid theoretical framework to study transitions, and a numerical procedure to characterize the most probable paths VDE1 ; VDE2 ; VDE3 ; kohn . The idea is to allow for fluctuactions that drive the system out of equilibrium near the transition threshold. The energy barriers between the minima are related to the times spent by the system in each basin of attraction, and to the probability of transition.

We prove that, under general hypotheses on the energy function, the conformational change occurs by a domino effect in which local destabilization events trigger their neighbors, and this propagates along the capsid until completion. This confirms the energy cascade hypothesis ross ; cardone and suggests that local interaction rules govern the details of the transition, which is henceforth not a collective, concerted motion of all capsomers, but rather propagates as a wave along the capsid.

The fact that the transition is an intrinsically local affair is also confirmed by linearizing the model and computing the stationary distribution of fluctuactions around metastable states: the concentration matrix shows that only adjacent pentamers are correlated. Explicit results for a special form of the capsid energy are presented and discussed.

II The dodecahedral model

According to the Caspar and Klug theory, actual viral capsids are icosahedral shells made of $T$-multiples of 60 proteins, with $T$ an integer CK . The proteins aggregate in small structural units called capsomers, that are usually composed by 2, 3, 5 or 6 proteins. The capsomers are almost rigid, and are often the basic subunits of assembly, as well as of disassembly. Motivated by bakker , we focus here on a class of viruses, such as Equine Rhynitis Virus A (ERAV), that release the genome by opening large pores in the capsid by translation and rotation of pentagonal rigid units made of 20 coat proteins. In this case, the capsid behaves as a dodecahedron made of 12 rigid pentagonal faces. The dodecahedon is also a popular model of viral capsid in theoretical investigations, since it is the smallest polyhedron exhibing icosahedral symmetry, still retaining a rich connectivity (each unit is coordinated to 5 other units).

Hence, in this paper we shall employ a dodecahedron $\mathcal{C}$ (see Figure 1) as a model of viral capsid. Consistently with this simplification, we assume that every configuration of the capsid is described by an order parameter ${\boldsymbol{x}}=(x_{1},\dots,x_{12})$, with indexing corresponding to a labeling of the faces as in Figure 1, where $x_{i}\in{\mathbb{R}}$ is a variable describing the state of each pentagonal unit. For instance, in destabilization problems leading to the opening of the capsid during the infection process, we can assume that each pentagonal face of the dodecahedron $\mathcal{C}$ is a rigid plaquette that can only translate along an axis orthogonal to its plane, and choose $x_{i}$ as its radial displacement.

In order to account for interactions between adjacent pentagons, it is convenient to work on the dual graph $G$ of $\mathcal{C}$, which has the property that the vertices of $G$ correspond to the faces of $\mathcal{C}$, the edges of $G$ to the edges of $\mathcal{C}$ and the faces of $G$ to the vertices of $\mathcal{C}$: $G$ is the graph of the icosahedron. From now on we view ${\boldsymbol{x}}$ as a field on the vertices of the icosahedral graph $G$.

Denoting by $V=\{1,\dots,12\}$ and $E$ the sets of vertices and edges of $G$, respectively, the adjacency matrix is the square symmetric $12\times 12$ matrix defined by

[TABLE]

Notice that the $ij$-entry of the adjacency matrix $A$ of $G$ is not vanishing if and only if the $i$ and $j$ faces of the dodecahedron $\mathcal{C}$ meet at a common edge.

II.1 Symmetry

We review here some basic relation between the adjacency matrix of $G$ and the symmetry of the associated polyhedron.

We say that a map $V\to V$ is an automorphism of $G$ if it is one-to-one and if it and its inverse maps adjacent vertices into adjacent vertices. The group of automorphisms of the icosahedral graph $G$ is the Coxeter group $\mathcal{H}_{3}={\mathcal{I}}\times\mathbb{Z}_{2}$ of order 120, where ${\mathcal{I}}$ denotes the rotational group of the icosahedron, with order 60 alg_graph . The group $\mathcal{I}$ acts on the vertices of $G$, which are the faces of $\mathcal{C}$, inducing a permutation representation (perm rep) $\sigma:{\mathcal{I}}\rightarrow S_{12}$, where $S_{12}$ is the symmetric group over 12 elements. The perm rep $\sigma$ induces a representation $\rho:{\mathcal{I}}\rightarrow GL(12,{\mathbb{R}})$, given by

[TABLE]

where ${\boldsymbol{e}}_{j},j=1,\ldots,12$ denotes the standard basis of ${\mathbb{R}}^{12}$.

The direct product decomposition of $\mathcal{H}_{3}$ implies that the representation

[TABLE]

is a representation of $\mathcal{H}_{3}$. Here $\Gamma=\{\pm 1\}$ is a representation of $\mathbb{Z}_{2}$ and $\otimes$ denotes the tensor product of representations (in this case, since the groups are finite, this is the Kronecker product of matrices) jones . Since $\mathcal{H}_{3}$ is the automorphism group of $G$, $A$ commutes with all the matrices of $\widetilde{\rho}(\mathcal{H}_{3})$.

There is a connection between the eigenspaces of $A$ and the decomposition into irreducible representations (irreps) of $\widetilde{\rho}$ spectra . The character table of ${\mathcal{I}}$ is given by

[TABLE]

where $\tau=\frac{1+\sqrt{5}}{2}$ is the golden ratio, $g_{2}$ and $g_{5}$ a two- and five-fold rotation of the icosahedron, respectively, and $C$ denotes the conjugacy class of an element of the group . The decomposition of $\widetilde{\rho}$ into irreps is given by ziac

[TABLE]

and there exists a matrix $R\in GL(12,{\mathbb{R}})$ such that $R^{-1}\widetilde{\rho}R=\widehat{\rho}$. The explicit form of $R$ is given in ziac . It is shown in spectra that the matrix $R$ diagonalises the adjacency matrix $A$ of the graph $G$; in particular, the spectrum of $A$ is given by

[TABLE]

III Energy

As discussed in the Introduction, we associate with the capsid an energy function $E:{\mathbb{R}}^{12}\rightarrow{\mathbb{R}}$ that is the sum of two contributions: a term driving the configurational change, that we identify with the internal pressure due to the confinement of the genomic material inside the capsid, and a term opposing the conformational change, that we identify here with the cohesive force between the pentamers at their interfaces. This assumption is consistent with the proposals of kegel2004 ; zandi_2 ; zandi2006 where it has been used to model assembly. Hence, we write the total energy of the capsid as

[TABLE]

with $f:{\mathbb{R}}\to{\mathbb{R}}$ and $g:{\mathbb{R}}^{2}\to{\mathbb{R}}$ to be specified below. The constant $\gamma$ is a real parameter controlling the stability of the capsid, and we assume that it depends on the environment of the capsid, such as pH, salt concentration, temperature.

Notice that, by construction, the energy is invariant under the action of the symmetry group of the capsid ziac . Writing, with a slight abuse of notation, $\widetilde{\rho}(\mathcal{H}_{3})\subset GL(12,{\mathbb{R}})$ simply as $\mathcal{H}_{3}$, this means that

[TABLE]

which implies that the gradient $\nabla E({\boldsymbol{x}})$ is equivariant, i.e.

[TABLE]

We require that $f:{\mathbb{R}}\to{\mathbb{R}}$ and $g:{\mathbb{R}}^{2}\to{\mathbb{R}}$ are smooth and

• $A_{1}$

the expansion energy is monotone decreasing and convex, and has an absolute minimum at state $x=1$, which we identify with the ’switched’ state of the pentamer (either fully detached of switched to the final conformation). Also, we assume that $f(1)=f^{\prime}(1^{-})=0$, $f^{\prime\prime}(1^{-})>0$. For $x>1$ the behavior of $f$ in this model is irrelevant: we can either assume that $f$ is constant, or that it is convex at $f=1$, which implies that the pentamer cannot detach further, or further modify its conformation. The states for which $x=1$ can be identified to the locked states of ross , i.e., intermediate states from which the transition cannot be reversed.

Notice that the convexity of $f$ as a function of a single variable does not imply the convexity of the energy $E$ as a function of the full state variable ${\boldsymbol{x}}$.

• $A_{2}$

the interaction energy $g=g(x,y)$ has a unique strict global mimimum at $x=y=0$, i.e., when both adjacent pentamers are attached to each other or in the initial state, and is symmetric, i.e., $g(x,y)=g(y,x)$. On the other hand, we must also require that the interaction becomes negligible when two pentamers are either sufficiently far, or one of them is in the switched state. The interaction radius is related to the intermolecular hydrophobic forces between capsomers, and should be much smaller than the distance at which Coulomb interaction is relevant kegel2004 .

Hence, writing $g_{1}=\frac{\partial g}{\partial x}$, $g_{2}=\frac{\partial g}{\partial y}$, $g_{11}=\frac{\partial^{2}g}{\partial x^{2}}$, and $g_{12}=\frac{\partial^{2}g}{\partial x\partial y}$, we require that $g_{1}(x,x)\geq 0$ and $g_{12}(x,x)\leq 0$ for every $x>0$, and there is a cutoff value $\bar{d}<<1$ such that $g$ vanishes identically for $x^{2}+y^{2}\geq\bar{d}^{2}$. In other words, $g(x,y)$ and all its derivatives vanish identically when one of the arguments equals 1. The hypothesis that $g_{12}(x,x)\leq 0$ implies that the transition pathway corresponding to $\gamma$ traversing the critical value (see next Section) is isotropic, consistent with normal mode analysis of capsid destabilization tama . It is possible to show that this requirement it is satisfied for the potential for hydrophobic interactions in kegel2004 .

We remark that both energies vanish identically when a component $x_{i}$ is equal to 1. As mentioned above, we identify states such that $x_{i}=1$ as states where the corresponding pentamer has switched to its final state, either by changing its conformation, or by detaching from the capsid, depending on the problem at hand. On the other hand, pentamers such that $x_{i}\sim 0$ will be viewed as being in the initial, or attached state, and will be said to be in state 0.

We also assume that every minimum of the energy is uniquely characterized by the corresponding combination of pentamers that are in state 1, i.e., by its combination of components that equal 1. Formally,

• $A_{3}$

if two minima of the energy $\hat{\boldsymbol{x}}$ and $\tilde{\boldsymbol{x}}$ are such that, for each $i$, either $\hat{x}_{i}=\tilde{x}_{i}=1$ or $\hat{x}_{i},\tilde{x}_{i}<1$, then $\hat{\boldsymbol{x}}=\tilde{\boldsymbol{x}}$.

• $A_{4}$

if two minima of the energy $\hat{\boldsymbol{x}}$ and $\tilde{\boldsymbol{x}}$ do not belong to the same icosahedral orbit, then $E(\hat{\boldsymbol{x}})\neq E(\tilde{\boldsymbol{x}})$: in other words, each minimum (modulo symmetry) is characterized by a unique energy level.

III.1 Icosahedral minima

In this section we partially characterize the changes of the energy landscape resulting from variations of the bond strength $\gamma$. We show that there is a critical threshold value for $\gamma$, such that below it the closed capsid is unstable. As $\gamma$, influenced by the environmental conditions, decreases below this critical value, the system undergoes a saddle-node bifurcation to an icosahedrally symmetric minimum which represents the final state. However, the analysis of the eigenvalues of the Hessian of the energy shows that the transition occurs by the activation of an icosahedrally symmetric ’breathing’ mode that corresponds to an icosahedral expansion of the capsid, a result inconsistent with the energy cascade hypothesis, but consistent with normal mode analysis of capsid destabilization tama .

First notice that assumption $A_{3}$ implies that there is at most equilibrium in which all pentamers are in the initial state, such that $x_{i}\neq 1$ for all $i$. Using the result (28) below, it follows that this equilibrium has full icosahedral symmetry, i.e., $x_{i}=x_{0}$, for all $i=1,\dots,12$.

The gradient and the Hessian of the energy are

[TABLE]

and

[TABLE]

Then $x_{0}$ must be a solution of the equation

[TABLE]

The above equation has always the solution $x=1$, i.e., the configuration in which all pentamers have switched to state 1. Further, by $A_{2}$, the function $g_{1}(x,x)$ is nonnegative and $g_{1}(0,0)=g_{1}(1,1)=0$. Since $f^{\prime}(x)<0$ for $x<1$, then there exists a critical value $\gamma_{c}$ such that for $\gamma>\gamma_{c}$ there are two solutions $x_{s}=x_{s}(\gamma)<1$ and $x_{u}=x_{u}(\gamma)<1$ of (7), with the property that $h^{\prime}(x_{s})>0$ and $h^{\prime}(x_{u})<0$, while for $\gamma<\gamma_{c}$ there are no solutions. Hence, the system undergoes a saddle-node bifurcation at $\gamma_{c}$.

The Hessian of the energy computed at an isotropic state $x_{i}=x$ for all $i$ is

[TABLE]

with $a=a(x)=f^{\prime\prime}(x)+10\gamma g_{11}(x,x)$ and $b=b(x)=2\gamma g_{12}(x,x)$. Notice that, by our hypotheses on $g,f$, we have that $b(x)\leq 0$, and the Hessian is proportional to the adjacency matrix $A$ of the graph (modulo the addition of a multiple of the identity). Therefore, its eigenspaces are also associated with the irreps of the representation $\widetilde{\rho}$ of $\mathcal{H}_{3}$ as in (1). The eigenvalues of the Hessian listed in increasing order are

[TABLE]

Hence,

[TABLE]

is the minimum eigenvalue of the Hessian, so that for $\gamma>\gamma_{c}$ the first critical point $x_{s}$ is a relative minimum of the energy, since the minimum eigenvalue $\mu_{1}(x_{s})$ of the Hessian is positive. A similar argument shows that $x_{u}$ is a saddle point. Notice also that, at the bifurcation point, the smallest eigenvalue of the Hessian changes sign, and the associated eigenspace is the isotropic line: hence, destabilization at the critical value $\gamma_{c}$ occurs by a breathing expansive mode.

IV Conformational changes

In many cases, the configurational changes of the capsid are triggered by variations of the chemical environment of the virion, which modifies the interactions between the proteins and, by consequence, the energy function. In our simplified model, such changes are reckoned by variations of the parameter $\gamma$.

However, the saddle-node bifurcation occurring at $\gamma_{c}$, in which the destabilization occurs via an icosahedrally symmetric expansion, does not account for the complexity of the interactions between the capsomers. In fact, the cooperative nature of the stability of the capsid suggests that the destabilization occurs as a cascade of elementary events: the switching of single capsomers requires less energy than the simultaneous change of state of all of them, and once one of them has switched, its neighbors are destabilized and a cascade of destabilization event is triggered, with decreasing energy barriers as the number of unswitched pentamers decreases.

A possible way to study the destabilization cascade is to work in the stable regime, and to treat the process in terms of a path visiting the metastable states of the energy according to the law of rare events, using Large Deviations Theory (LDT) for stochastic dynamical systems. In other terms, we assume that, near the initial-state equilibrium, additive fluctuactions of the system due to changes of the environment are able to drive it out of equilibrium, with increasing probability as $\gamma\to\gamma_{c}$, and the actual conformational change can be described by this stochastic process.

Hence, we assume from now on that $\gamma>\gamma_{c}$, and denote by ${\boldsymbol{x}}_{0}=(x_{s},\dots,x_{s})$ the icosahedral minimum of the energy $E$ for which $x_{s}<1$, and by ${\boldsymbol{x}}_{1}=(1,\dots,1)$ the final-state icosahedral minimum. In problems involving the opening of the capsid, ${\boldsymbol{x}}_{0}$ and ${\boldsymbol{x}}_{1}$ correspond to the closed and open capsid, respectively.

One way to study the stability of a complex structure is to construct a dynamical system whose attractors are the local minima of the energy, i.e., to assume that deviations from equilibria satisfy a deterministic dissipative gradient-flow dynamics

[TABLE]

where the superposed dot denotes differentiation with respect to time $t\in{\mathbb{R}}^{+}$, and ${\boldsymbol{x}}={\boldsymbol{x}}(t)$ is the motion in the configuration space ${\mathbb{R}}^{12}$. Since $\gamma>\gamma_{c}$, the system (8) has an attractor at ${\boldsymbol{x}}={\boldsymbol{x}}_{0}$ (i.e., $x_{i}=x_{0}=x_{s}$ $\forall i$): the initial state is stable for the dynamics (8), and every small perturbation of this state tends to vanish over finite time intervals.

Now, let ${\boldsymbol{W}}=(W^{(1)}(t),\ldots,W^{(12)}(t))$ be a $12$-dimensional Brownian motion defined on a probability space $(\Omega,\mathcal{A},\mathbb{P})$. We consider random perturbations of the dynamical system (8), by superimposing an additive noise $dW^{(i)}$ to each face $i$ of the polyhedron ${\mathcal{C}}$. This results in the stochastic differential equation

[TABLE]

with $\epsilon>0$ a small parameter. We remark that the Langevin dynamics (9) is here regarded only as a tool to explore the energy landscape of the system: the actual dynamics of the conformational change must involve more refined physical considerations. In the overdamped Langevin equation $\epsilon^{2}$ is proportional to the temperature.

IV.1 Results from large deviations theory

In this section we briefly review some key concepts from LDT that we are going to use throughout the paper. We refer to freidlin for the general theory.

The main object of LDT is the action functional $S_{T}({\boldsymbol{\varphi}})$ which, for the equation (9), is given by

[TABLE]

where $T>0$ and ${\boldsymbol{\varphi}}:[0,T]\rightarrow{\mathbb{R}}^{12}$ is an absolutely continuous path (actually, the functional is defined on continuous paths, but its value is $+\infty$ if ${\boldsymbol{\varphi}}$ is not absolutely continuous). The LDT states that, given a bounded set $D\subseteq{\mathbb{R}}^{12}$, then the probability ${\mathbb{P}}^{{\boldsymbol{x}}}({\boldsymbol{x}}(T)\in D)$ that the solution ${\boldsymbol{x}}(T)$ of (9) belongs to $D$ given that ${\boldsymbol{x}}(0)={\boldsymbol{x}}$ satisfies

[TABLE]

where ${\mathcal{C}}_{{\boldsymbol{x}},D}=\{{\boldsymbol{\varphi}}\in{\mathcal{C}}([0,T],{\mathbb{R}}^{12}):{\boldsymbol{\varphi}}(0)={\boldsymbol{x}},{\boldsymbol{\varphi}}(T)\in D\}$. If the event occurs, then ${\boldsymbol{x}}(t)$ is arbitrarily close to the minimizer

[TABLE]

in the sense that, for every $\delta>0$,

[TABLE]

Another central object in LDT is the quasipotential

[TABLE]

where ${\mathcal{C}}_{{\boldsymbol{x}},{\boldsymbol{y}}}=\{{\boldsymbol{\varphi}}\in{\mathcal{C}}([0,T],{\mathbb{R}}^{12}):{\boldsymbol{\varphi}}(0)={\boldsymbol{x}},{\boldsymbol{\varphi}}(T)={\boldsymbol{y}}\}$. For large time intervals, if ${\boldsymbol{x}}_{0}$ is the unique stable point of (8), the density $p({\boldsymbol{x}})$ associated with the stationary distribution of (9), provided it exists, is such that

[TABLE]

where $Z$ is a normalization constant. Here $\asymp$ denotes the log-asymptotic equivalence for $\epsilon\to 0$, i.e. $\lim_{\epsilon\to 0}\epsilon^{2}\ln(p({\boldsymbol{x}}))=-V({\boldsymbol{x}}_{0},{\boldsymbol{x}})$. For a gradient system like (9), the quasipotential at a point ${\boldsymbol{x}}$ lying in the basin of attraction of a minimum ${\boldsymbol{x}}_{0}$ of $E$ is given by

[TABLE]

Suppose ${\boldsymbol{x}}_{A}$ and ${\boldsymbol{x}}_{B}$ are two (local) minima of $E$ separated by a single saddle point ${\boldsymbol{x}}_{S}$. Then the minimizer ${\boldsymbol{\varphi}}^{*}$ of the action functional $S_{T}({\boldsymbol{\varphi}})$ with ${\boldsymbol{\varphi}}(0)={\boldsymbol{x}}_{A}$ and ${\boldsymbol{\varphi}}(T)={\boldsymbol{x}}_{B}$ is the path such that, up to a normalization constant,

[TABLE]

where $(\nabla E)^{\perp}$ is the component of $\nabla E$ normal to ${\boldsymbol{\varphi}}^{*}$:

[TABLE]

where $\widehat{{\boldsymbol{\tau}}}$ is the unit tangent of ${\boldsymbol{\varphi}}^{*}$ and $\langle,\rangle$ denotes the standard Euclidean product. The path ${\boldsymbol{\varphi}}^{*}$ is referred to as the Minimum Energy Path (MEP) connecting ${\boldsymbol{x}}_{A}$ and ${\boldsymbol{x}}_{B}$ VDE1 .

The dynamics of (9) on the energy landscape can be modeled as a discrete time Markov chain with states the minima of $E({\boldsymbol{x}})$ freidlin ; cameron1 . In this framework, two minima ${\boldsymbol{x}}_{\alpha}$ and ${\boldsymbol{x}}_{\beta}$ are connected if the minimizer of $S_{T}({\boldsymbol{\varphi}})$ corresponds to the MEP ${\boldsymbol{\varphi}}^{*}$ with a single maximal value of the energy along it. The off-diagonal entries of the generator matrix $Q$ of the chain are given by

[TABLE]

where ${\boldsymbol{x}}_{\alpha\beta}$ is the unique saddle point between ${\boldsymbol{x}}_{\alpha}$ and ${\boldsymbol{x}}_{\beta}$. The diagonal entries are chosen such that the sum of each row is zero:

[TABLE]

The jump matrix $\Pi$ is obtained from $Q$ by setting

[TABLE]

The jump matrix $\Pi$ is the generator matrix of a discrete Markov chain prob . Following cameron1 , we introduce the limiting jump matrix (or zero-temperature jump matrix)

[TABLE]

The limiting jump matrix depends only on the values of the potential at the saddles cameron1 .

IV.2 Admissible states

The energy landscape is complex and there are many local minima, among which there is the icosahedral minimum. By assumption $A_{3}$, local minima are completely characterized by the combination of pentamers $i$ such that $x_{i}=1$, and $0<x_{j}<<1$ for the remaining indices, i.e., by the combination of detached pentamers. In this section we make precise this notion.

We parametrize the minima by strings ${\boldsymbol{s}}\in\{0,1\}^{12}$, where $s_{i}=1$ means that the pentagon $i$ is in state 1, for instance detached from the capsid, while $s_{i}=0$ means that the pentagon is still in its initial state, for instance attached to its neighbors. Precisely, we say that a state ${\boldsymbol{s}}$ is admissible if it belongs to the basin of attraction of a critical point $\hat{\boldsymbol{x}}({\boldsymbol{s}})$ of the dynamical system (8) that corresponds to the same combination of pentamers in state 1 as ${\boldsymbol{s}}$, i.e., if

[TABLE]

where ${\boldsymbol{x}}(t,{\boldsymbol{s}})$ is the solution of (8) with initial conditions ${\boldsymbol{x}}(0)={\boldsymbol{s}}$.

We denote by $S$ the set of admissible states, and write $\sum_{i}s_{i}$ for the number of pentamers in state 1. Special admissible states are the initial configuration ${\boldsymbol{s}}_{0}=(0,0,\dots,0)$ and the final configuration ${\boldsymbol{s}}_{1}=(1,1,\dots,1)$.

As an example, a set of admissible states for the energy (39) is shown in Figure 4. These are obtained by solving numerically the dynamical system (8) with initial conditions ${\boldsymbol{x}}(0)={\boldsymbol{s}}$, to determine the attractor $\hat{x}({\boldsymbol{s}})$ whose basin ${\boldsymbol{x}}(0)$ belongs to.

As seen in Section IV.1, the dynamics on the energy landscape can be analyzed by defining a Markov chain on the $S$ of admissible states. This analysis is carried out in the next section.

IV.3 Transitions

In this section we introduce the basic tool that will allow to study transitions, i.e., the notion of connectivity between admissible states. The first basic idea here is that two states are connected if it is possible to go from one to the other by detaching pentamers. This automatically forbids reversible transitions in which pentamers can reattach or switch back to their initial state. This hypothesis is consistent with the notion of ’locked’ states proposed for the maturation pathway of HK97 ross , but is also reasonable in the disgregation problem since detached capsomers cannot usually reattach to the capsid.

The second basic notion is that two states are connected if the minimum energy path between them does not visit other minima, which is also a reasonable assumption here.

Precisely, given two admissible states ${\boldsymbol{s}},{\boldsymbol{s}}^{\prime}\in S$, we say that ${\boldsymbol{s}}<{\boldsymbol{s}}^{\prime}$ if $s_{i}^{\prime}=1$ when $s_{i}=1$ and $\sum_{i}s_{i}<\sum_{i}s_{i}^{\prime}$: in this case, since pentamers that are in state 1 in ${\boldsymbol{s}}$ are also in state 1 in ${\boldsymbol{s}}^{\prime}$, the state ${\boldsymbol{s}}^{\prime}$ is obtained by switching to $1$ some of the pentamers that are in state [math] in ${\boldsymbol{s}}$. The relation $<$ induces a strict partial order on the set $S$ consistent with the structure of directed graph below.

We say that there is a directed edge between two admissible states ${\boldsymbol{s}}$, and ${\boldsymbol{s}}^{\prime}$, and write ${\boldsymbol{s}}\to{\boldsymbol{s}}^{\prime}$, if

• i)

${\boldsymbol{s}}^{\prime}>{\boldsymbol{s}}$, so that no reverse transition or reattachment of pentamers is allowed; and

• ii)

there is a MEP connecting the corresponding minima $\hat{\boldsymbol{x}}({\boldsymbol{s}}),\hat{\boldsymbol{x}}({\boldsymbol{s}}^{\prime})$, along which the energy has a single maximum, as described in Section IV.1.

We remark again that requirement (i) is strongly restrictive, in that it excludes conformational changes involving pentamers going back to their initial configuration, or reattaching to the capsid. According to LDT, such transitions are indeed possible in the stochastic dynamics (9) and this may give rise to cycles. The Markov chain that we construct below, however, is meant to describe a restricted situation in which configurational changes are irreversible, as is the case for most configurational changes in capsids, first of all the disgregation. In other terms, the values $x_{i}=1$ act as absorbing states for the $i$-th component of SDE (9).

If ${\boldsymbol{x}}_{m}$ is the point of the MEP between $\hat{\boldsymbol{x}}({\boldsymbol{s}})$ and $\hat{\boldsymbol{x}}({\boldsymbol{s}}^{\prime})$ where the maximum is attained, we define the barrier between ${\boldsymbol{s}}$ and ${\boldsymbol{s}}^{\prime}$ by

[TABLE]

which is positive by construction. We set $\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})=+\infty$ if ${\boldsymbol{s}}$ is not connected to ${\boldsymbol{s}}^{\prime}$.

The above procedure endows the set $S$ with a structure of directed acyclic graph with positive weights. To construct explicitly the weights $\beta$, the MEP between any pair of states can be determined using the numerical procedure introduced in VDE1 ; VDE2 ; VDE3 , taking for instance as initial path the straight line joining $\hat{\boldsymbol{x}}({\boldsymbol{s}})$ to $\hat{\boldsymbol{x}}({\boldsymbol{s}}^{\prime})$.

IV.4 The Markov chain

In order to describe the dynamics on the energy landscape and determine the most probable transitions paths, as discussed in Section IV.1 we construct a Markov chain on the state space $S$. The generator matrix $Q$ of this chain is given by (compare with (17))

[TABLE]

and

[TABLE]

The jump matrix $\Pi$ as in (18) is then given by

[TABLE]

Letting $\epsilon\to 0$ we obtain the zero-temperature jump matrix (cf. (19))

[TABLE]

The matrix $\Pi_{0}$ defines a Markov chain on the set of admissible states $S$ by assuming that the transition probability between the states ${\boldsymbol{s}},{\boldsymbol{s}}^{\prime}$ is

[TABLE]

The explicit form of the matrix $\Pi_{0}$ is the following:

[TABLE]

where $N({\boldsymbol{s}})=|\{{\boldsymbol{s}}^{\prime}:\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})=\min_{{\boldsymbol{s}}\to{\boldsymbol{s}}^{\prime\prime}}\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime\prime})\}|$ is the number of states that can be reached from ${\boldsymbol{s}}$ along a path with minimum barrier.

The most probable transition path between the closed state ${\boldsymbol{s}}_{0}$ to the open state ${\boldsymbol{s}}_{1}$ is the realization of the chain that maximizes the transition probability at each step, i.e., the path along which the barriers are minimal among all the admissible transitions outgoing from each vertex (cf. Theorem 6.6.1 in freidlin ).

We point out that the use of the zero-temperature matrix $\Pi_{0}$ in the construction of the Markov chain implies that we are only considering transitions between nearest neighboring states (which correspond to the first order cycles in the terminology of freidlin ; cameron1 ). The construction of higher order cycles would allow a wider analysis of the transitions between pentamers. However, the configuration space for the latter would become extremely complex, and the analysis carried out here would be unfeasible.

V Reduction by symmetry

In order to reduce the complexity of the problem, we will now restate the above results in terms of symmetry classes of minima. In fact, since the energy $E$ is invariant under the action of the symmetry group $\mathcal{H}_{3}$, minima are mapped into minima by $\mathcal{H}_{3}$, and this induces a permutation action on the set of minima of $E$. By (4) and uniqueness of the solution of (8), if ${\boldsymbol{x}}(t)$ is the solution of (8) with initial condition ${\boldsymbol{s}}$, then $H{\boldsymbol{x}}(t)$ is the solution of (8) with initial condition $H{\boldsymbol{s}}$, for every $H\in\mathcal{H}_{3}$, and hence, by (20), it follows that

[TABLE]

so that, in turn, $\mathcal{H}_{3}$ acts on $S$. Denoting by

[TABLE]

the isotropy group of ${\boldsymbol{s}}$, (26) implies that $\hat{\boldsymbol{x}}({\boldsymbol{s}}))$ is invariant under $\text{Fix}({\boldsymbol{s}})$, i.e.,

[TABLE]

which in turn implies that $\hat{\boldsymbol{x}}({\boldsymbol{s}}_{0})$ and $\hat{\boldsymbol{x}}({\boldsymbol{s}}_{1})$ have all components equal, since $\text{Fix}({\boldsymbol{s}}_{0})=\text{Fix}({\boldsymbol{s}}_{1})=\mathcal{H}_{3}$ and therefore

[TABLE]

The set $S$ therefore can be decomposed into orbits of $\mathcal{H}_{3}$. We denote by $\Sigma$ the set $S/\mathcal{H}_{3}$. Figure 4 shows the reduced state space $\Sigma$ for the energy (39). We now show that the Markov chain $\Pi$ on $S$ induces a Markov chain $\tilde{\Pi}$ on $\Sigma$.

Consider first a MEP ${\boldsymbol{\varphi}}^{*}$ connecting two minima ${\boldsymbol{x}}_{0}$ and ${\boldsymbol{x}}_{1}$: by construction, ${\boldsymbol{\varphi}}^{*}$ is a minimizer of the action functional $S_{T}({\boldsymbol{\varphi}})$ as in (10). Since $E$ is invariant under $\mathcal{H}_{3}$, by (4) $\nabla E(H{\boldsymbol{x}})=H\nabla E({\boldsymbol{x}})$ for $H\in\mathcal{H}_{3}$, and therefore, since $H$ is orthogonal, $S_{T}({\boldsymbol{\varphi}})=S_{T}(H{\boldsymbol{\varphi}})$. This means that if ${\boldsymbol{\varphi}}^{*}$ is a MEP from ${\boldsymbol{x}}_{0}$ to ${\boldsymbol{x}}_{1}$, then $H{\boldsymbol{\varphi}}^{*}$ is a MEP from $H{\boldsymbol{x}}_{0}$ to $H{\boldsymbol{x}}_{1}$. Hence, noting that ${\boldsymbol{s}}<{\boldsymbol{s}}^{\prime}$ implies $H{\boldsymbol{s}}<H{\boldsymbol{s}}^{\prime}$,

[TABLE]

We now prove a basic property of barriers. From the invariance of $E$ and the above discussion, it follows that $E({\boldsymbol{\varphi}}(s))=E(H{\boldsymbol{\varphi}}(s))$ for every $s\in[0,T]$, and using (26) and (21) this shows that

[TABLE]

Notice in particular that, denoting by $\text{Fix}({\boldsymbol{s}})=\left\{H\in\mathcal{H}_{3}:H{\boldsymbol{s}}={\boldsymbol{s}}\right\}$ the isotropy group of ${\boldsymbol{s}}$,

[TABLE]

Hence, the value of the barrier between the state ${\boldsymbol{s}}$ and all states in $\text{Fix}({\boldsymbol{s}}){\boldsymbol{s}}^{\prime}$ is the same.

We define a weighted graph with vertex set $\Sigma$ as follows: for ${\boldsymbol{\sigma}},{\boldsymbol{\sigma}}^{\prime}\in\Sigma$, we write

[TABLE]

if there exist ${\boldsymbol{s}}\in{\boldsymbol{\sigma}}$, ${\boldsymbol{s}}^{\prime}\in{\boldsymbol{\sigma}}^{\prime}$ such that

[TABLE]

in which case we write

[TABLE]

The above definition is well given. In fact, let ${\boldsymbol{s}}^{\prime}\in{\boldsymbol{\sigma}}^{\prime}$ such that $\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})=\min_{{\boldsymbol{s}}^{\prime\prime}\in{\boldsymbol{\sigma}}^{\prime}}\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime\prime})$, and assume that there exist ${\boldsymbol{t}}\in{\boldsymbol{\sigma}}$, ${\boldsymbol{t}}^{\prime}\in{\boldsymbol{\sigma}}^{\prime}$ such that $\beta({\boldsymbol{t}},{\boldsymbol{t}}^{\prime})=\min_{{\boldsymbol{t}}^{\prime\prime}\in{\boldsymbol{\sigma}}^{\prime}}\beta({\boldsymbol{t}},{\boldsymbol{t}}^{\prime\prime})\neq\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})$. Then by transitivity on the orbits there exists $H\in\mathcal{H}_{3}$ such that ${\boldsymbol{t}}=H{\boldsymbol{s}}$, so that by (30) $\beta({\boldsymbol{t}},{\boldsymbol{t}}^{\prime\prime})=\beta({\boldsymbol{s}},H^{\top}{\boldsymbol{t}}^{\prime})$. Clearly, $H^{\top}{\boldsymbol{t}}^{\prime}\in{\boldsymbol{\sigma}}^{\prime}$ and this implies that $\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})\leq\beta({\boldsymbol{t}},{\boldsymbol{t}}^{\prime\prime})$. Since the converse also holds, we arrive at a contradiction and $\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})=\beta({\boldsymbol{t}},{\boldsymbol{t}}^{\prime})$.

The above construction defines a Markov chain $\tilde{\Pi}$ on $\Sigma$:

[TABLE]

To compute the low-temperature limit of (35), we make a further assumption on the system:

• $A_{5}$

For every admissible state ${\boldsymbol{s}}$, if $\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime})=\beta({\boldsymbol{s}},{\boldsymbol{s}}^{\prime\prime})$, there exists $H\in\mathcal{H}_{3}$ such that ${\boldsymbol{s}}^{\prime\prime}=H{\boldsymbol{s}}^{\prime}$.

This hypothesis is consistent with assumption $A_{4}$ on the minima of the energy, and is necessary to further reduce the complexity of the problem, since it guarantees that, starting from a given state, there is a unique orbit to which the system can transform along a MEP.

Now, as $\epsilon\to 0$, the dominant terms at the denominator of (35) correspond to those ${\boldsymbol{\sigma}}^{\prime\prime}$ such that $\tilde{\beta}({\boldsymbol{\sigma}},{\boldsymbol{\sigma}}^{\prime\prime})=\min_{{\boldsymbol{\sigma}}\to{\boldsymbol{\sigma}}^{\prime\prime\prime}}\tilde{\beta}({\boldsymbol{\sigma}},{\boldsymbol{\sigma}}^{\prime\prime\prime})$. By assumption $A_{5}$, there can be only one orbit that realizes the minimum, so that the leading term only contains one summand. Hence, the low-temperature limit of (35) is

[TABLE]

VI Analysis of a special model

We present below the analysis of a specific model, meant to describe the opening and disgregation of the capsid from the initial closed state and a final state in which all pentamers have detached from the capsid, motivated by the putative mechanism by which ERAV releases its genome within the host cell bakker . In this context, $x_{i}$ is the radial displacement of the capsomers along their axes (cf. Fig. 2). The model is based on the special form of energy function

[TABLE]

where $f$ is the expansion energy and $g$ is the attractive interaction energy, that depends on the squared distance $d:{\mathbb{R}}^{2}\rightarrow{\mathbb{R}}$ between the capsomers, defined by

[TABLE]

where $\alpha$ is the angle between two neighboring icosahedral axes (cf. Figure 2).

Precisely, we choose the expansion and attractive energies, and the constant $\gamma$ as

[TABLE]

The expansion energy $f$ vanishes for $x\geq 1$, while the attractive energy $g$ has a much shorter radius. The attractive energy with fast exponential increase is consistent with the form proposed in kegel2004 (formula 1) for hydrophobic interactions between apolar surfaces, while the expansive energy is just the simplest choice of a polynomial monotone decreasing convex function consistent with hypothesis $A_{1}$.

For each given number of detached pentamer, a representative of each icosahedral orbit of admissible states is shown in Fig. 4. Inspection of that figure shows that

• •

There are multiple states with decreasing number of attached pentagons, but there is only one state (modulo the icosahedral group) with 5 and 4 attached pentamers.

• •

There are no admissible states with less than 4 attached pentamers. While it is obvious that a single pentamer cannot be attached to anything, even clusters of 4 or more pentamers cannot be stable in real capsids. This is a consequence of working with such a simplified model.

• •

In each admissible state the pentamer configuration is connected.

In Tables 1 and 2 below are listed the half-barriers $\beta({\boldsymbol{\sigma}},{\boldsymbol{\sigma}}^{\prime})/2$ between symmetry classes of admissible states for the detachment problem and energy (39). Only the half-barriers between nearest-neighbor states, that correspond to the detachment of a single pentamer, are indicated, with the exception of some barriers between states that differ by at least two pentamers. These are written only when they are lower than the barriers between nearest neighbors, or when there is no stable nearest neighbor.

Inspection of Tables 1 and 2 shows that the lowest barrier between a state and its out neighbors mostly corresponds to the detachment of one of the pentamers with the lowest number of attached neighbors. Computing the non-zero entries of the zero-temperature shows that the most probable transition path from the closed to the open state is

[TABLE]

withe the same labels as in Figure 4. Figure 3 shows a realization of the most probable transition path.

VII Fluctuactions near equilibrium

Linearization of the dynamical system (9) near equilibria gives information on the correlation between the fluctuactions of the building blocks of the system.

Consider an admissible state ${\boldsymbol{s}}$ and the corresponding equilibrium $\hat{\boldsymbol{x}}({\boldsymbol{s}})$. The Hessian

[TABLE]

of the energy at this point is positive definite by assumption (cf. (20)), and, using (3), (4) and (26), we see also that

[TABLE]

Fluctuations near the stable equilibrium $\hat{\boldsymbol{x}}({\boldsymbol{s}})$ are ruled by the linearized system

[TABLE]

Assuming the initial condition to be constant or normally distributed, the explicit solution of (42) is a Gaussian process given by (see arnold )

[TABLE]

The expectation ${\boldsymbol{m}}={\boldsymbol{m}}(t):=\mathbb{E}[{\boldsymbol{x}}(t)]$ and the covariance matrix $P(t)$ of ${\boldsymbol{x}}(t)$, are solutions of the linear ODEs

[TABLE]

The solution of (44)1 is immediate and it is given by

[TABLE]

We find the explicit solutions of (44)2 using some results from matrix algebra. We associate with a matrix $B\in GL(n,{\mathbb{R}})$ the (column) vector $v(B)\in{\mathbb{R}}^{n^{2}}$ of its entries:

[TABLE]

The following property holds for every $B,C,D\in GL(n,{\mathbb{R}})$ barnett :

[TABLE]

where $\otimes$ denotes the Kronecker product of matrices. We transform (44)2 into a vector equation obtaining

[TABLE]

which is a linear non homogenous first order ODE, whose solutions is given by, setting $\widetilde{K}=K\otimes I+I\otimes K$,

[TABLE]

Since $K\otimes I$ and $I\otimes K$ commute, we have, using the properties of the Kronecker product (see johnson ),

[TABLE]

Using properties (46), equation (47) becomes

[TABLE]

Hence, the solution of (44)2 is given by

[TABLE]

The solution ${\boldsymbol{x}}(t)$ of (42) has therefore distribution $\mathcal{N}\left({\boldsymbol{m}}(t),P(t)\right)$. Since $K$ is positive definite, we have

[TABLE]

Hence, ${\boldsymbol{x}}(t)$ converges in distribution to a Gaussian random variable ${\boldsymbol{X}}$ with mean ${\boldsymbol{m}}={\boldsymbol{0}}$ and covariance matrix $Q=\frac{\epsilon^{2}}{2}K^{-1}$.

We point out that this result agrees with the LDT framework. In fact, the quasipotential of (42) as in (12) in the basin of attraction of ${\boldsymbol{x}}_{0}$ is given by $V({\boldsymbol{x}})={\boldsymbol{x}}\cdot K{\boldsymbol{x}}$. It follows from (13) that the stationary distribution $p({\boldsymbol{x}})$ of (42) is asymptotically given by

[TABLE]

which is the density of a Gaussian distribution $\mathcal{N}\left(\mathbf{0},\frac{1}{2}\epsilon^{2}K^{-1}\right)$. Here $Z$ is a normalization constant.

The limiting distribution $p({\boldsymbol{x}})$ has important consequences. Recall that, given three continous random variables $X$, $Y$ and $Z$ with joint density distibution $p_{XYZ}(x,y,z)$, we say that $X$ is conditionally independent of $Y$ given $Z$, and write $X\perp Y|Z$, if and only if

[TABLE]

Let ${\boldsymbol{X}}=(X_{1},\ldots,X_{n})\sim\mathcal{N}({\boldsymbol{m}},Q)$ be a multivariate normal distribution. The concentration matrix $C$ of ${\boldsymbol{X}}$ is the inverse of the covariance matrix $Q$ (provided that $\det(Q)\neq 0$). The entries of $C$ measure the correlation between the components. In particular, the partial correlation coefficients are given by

[TABLE]

where $\mathcal{S}=\{1,\ldots,n\}$. Moreover, the following holds:

[TABLE]

If condition (51) is satisfied, then ${\boldsymbol{X}}$ is known as a Gaussian Markov random field graphmodels . In our case, the concentration matrix of ${\boldsymbol{x}}$ is $\frac{2}{\epsilon^{2}}K$, and we have, from (50)

[TABLE]

Now, from the expression (6), it follows that

[TABLE]

since, by $A_{2}$, $g(x,y)$ and all its derivatives vanish whenever one of its arguments is 1. Moreover, if $\hat{x}_{i}<1$, so that $s_{i}=0$, then

[TABLE]

A first consequence of (6), (52) and (53) is that, at an admissible state, two pentamers are correlated if and only if they are adjacent in the configuration ${\boldsymbol{s}}$. In other words, the random variable ${\boldsymbol{x}}$ is a Gaussian graphical model with graph the subgraph of the icosahedral graph induced by the vertices with zero components of ${\boldsymbol{s}}$.

In turn, (54) shows that the diagonal elements $K_{ii}$ of the Hessian depend only on the pentamers adjacent to $i$ in the state ${\boldsymbol{s}}$. Actually, (54) suggests that $K_{ii}$ is larger the greater is the connectivity of pentamer $i$ in the state ${\boldsymbol{s}}$. Hence, a pentamer with high connectivity has a small correlation coefficient with its neighbors.

VIII Conclusions

Large deviations theory for stochastic differential equations is based on the notion that arbitrarily small stochastic perturbations are able to lead any system out of equilibrium, over sufficiently long times, and minimum energy paths allow to determine the most probable transition paths between metastable states: the resulting Markov chain on the set of minima completely describes the stochastic dynamics on the energy landscape. In this work we have used this approach to get insights into the process by which viral capsids change configuration, in either maturation or infection.

Our analysis supports the conjecture that destabilization occurs by a cascade of local events cardone ; ross : in fact, the energy is the sum of a destabilizing term on each pentamers, which does not depend on its connectivity, and a cohesive term opposing the transition, which accounts for the interactions between adjacent pentamers. Destabilization is the result of the competition between these terms: the expansive energy decreases whenever a pentamer changes state, but this requires providing an amount of cohesive energy proportional to the number of bonds broken in the process. This suggests that the cascade occurs by destabilization of those pentamer that have less bonds, and therefore have to pay less energy to detach. Simulations confirm the intuitive picture above: for instance, in the parameter range in which we are working, whenever a pentamer has only one bond left, it switches to its final configuration, because the energy gained by this process is larger than the amount lost by bond breaking.

For simplicity, we have computed in this paper just the zero-temperature transitions, i.e., the limit of the transition matrix of the Markov chain as $\epsilon\to 0$, but the procedure is fully general and allows one to study all transitions between metastable states.

In order to further explore how random fluctuactions affect the stability of a complex interacting structure, we have also focused on the linearization of the system near a metastable state. Clearly, the destabilization dynamics is strongly nonlinear, but the analysis of its linearization in the neighborhood of an attractor yields interesting information. First of all, the process is Gaussian and its limit distribution is a multivariate normal. The concentration matrix $C$ is proportional to the Hessian of the energy at equilibrium, which has nonzero entries only when two pentamers are adjacent. This is a first confirmation of the effect of locality on this model, but the point is that this allows to compute the conditional correlation coefficients between fluctuactions at adjacent pentamers. This shows that if a pentamer $i$ has many unbroken bonds and $H_{ii}$ and $C_{ii}$ are large, the fluctuactions are highly concentrated at that pentamer, and the correlation between this and adjacent pentamers is small.

On the other hand, when a pentamer $i$ has a small connectivity, the concentration coefficient $C_{ii}$ is small, and the pentamer has large negative correlation coefficients with its neighbors. This means that fluctuactions tend to amplify, which is a clue of the destabilization effect.

The procedure described in this paper can be automatized and generalized to any system made of pairwise interacting building blocks, once the interactions are encoded in a simple graph.

Acknowledgements.

We thank R. Kohn for suggesting the use of LDT to study transitions, and E. Vanden-Eijnden for valuable comments and suggestions. GI acknowledges support by the Italian GNFM (Progetto Giovani 2016) and GI and PC acknowledge support by the by the University of Torino (research project ’Modelli Aleatori’)).