Statistical Mechanical Treatments of Protein Amyloid Formation

Protein aggregation is an important field of investigation because it is closely related to the problem of neurodegenerative diseases, to the development of biomaterials, and to the growth of cellular structures such as cyto-skeleton. Self-aggregation of protein amyloids, for example, is a complicated process involving many species and levels of structures. This complexity, however, can be dealt with using statistical mechanical tools, such as free energies, partition functions, and transfer matrices. In this article, we review general strategies for studying protein aggregation using statistical mechanical approaches and show that canonical and grand canonical ensembles can be used in such approaches. The grand canonical approach is particularly convenient since competing pathways of assembly and dis-assembly can be considered simultaneously. Another advantage of using statistical mechanics is that numerically exact solutions can be obtained for all of the thermodynamic properties of fibrils, such as the amount of fibrils formed, as a function of initial protein concentration. Furthermore, statistical mechanics models can be used to fit experimental data when they are available for comparison.

perspectives, it is not surprising that the papers on protein aggregation differ widely in their emphasis and methodologies: from fundamental research related to molecular mechanisms and aggregation pathways to searching for biomarkers, drug targets, even to imaging of plaques in the brain, dissolution of fibrils and amyloids in vivo, etc. The present article is only concerned with the fundamental investigations into the aggregation mechanisms of amyloid formation related to neurodegenerative disease. Many review articles have been written on the approaches based on molecular dynamical simulations [1][2][3][4], as well as kinetic studies [5][6][7]. In our examination, we will focus instead on the statistical mechanical approaches to equilibrium assembly processes and present some new results while summarizing past approaches.
Early applications of statistical mechanical methods to the studies of protein problems can best be represented by the treatment of helix-coil transitions in proteins by Zimm and Bragg in the 1950s [8][9][10][11][12][13]. They assumed that each peptide bond linking amino acid residues together can exist in two states: a helical or non-helical state, and characterized the linear chain of residues with a partition function, which is a sum of all possible combinations of states.
It is in a similar spirit that our statistical mechanical treatment of protein aggregation has been developed [21], which is the main subject of this article. Our formalism of the aggregation processes was stimulated by other statistical mechanical studies. These works will be briefly reviewed in Section 3, along with conceptual developments of statistical mechanical techniques beyond those used in the Zimm-Bragg model. In Sections 4-6, canonical ensemble and grand canonical ensemble treatments of the aggregation processes will be separately presented, which is followed by a conclusion section. In Section 2, we first review some properties of amyloid proteins and aggregates.

Amyloid Aggregation
To see the complication of protein aggregation processes, we use β-amyloid as an example. Under proper conditions, such as higher concentration, amyloid monomers can aggregate into dimers, trimers, tetramers, . . . , oligomers. These co-existing oligomers are in rapid kinetic equilibrium, making it difficult to determine their structures or numbers [27]. Oligomers of the same size may exist in different conformations: some are partially ordered and some are disordered. As soluble oligomers grow larger they can become richer in β-sheet structures, but overall they lack secondary structures. They may resemble micelles, which have a spherical or cylindrical shape [28]. Aβ oligomers seem to range in size, with a diameter ranging from 5-15 nm and molar mass ranging from 20-50 kDa up to 1 MDa [29,30]. Additionally, oligomers with ordered β-sheet or β-hairpin structures are believed to form protofibrils at a higher rate than their disordered counterparts [3,[31][32][33]. For example, Li et al. [3] showed in numerical simulations that a native chain has to unfold partially into an intermediate with beta-hairpin structure before ordered assembly can be formed. Since protein aggregation is thought to be a nucleation process, some of the ordered oligomers may act as paranuclei [34]. Once formed, paranuclei can lead to the formation of protofibrils in down-hill fashion. The nucleation is illustrated in Figure 1b,c. Since Aβ oligomers may contain β-sheet structure, they could be pathway intermediate [28]. However, the same cannot be said about other oligomers, such as those comprised of prion proteins [35]. Additionally, Aβ  and Aβ  monomers may self-associate to form off-pathway globular assemblies including amylospheroids [34,36] and β-amyloid balls [37] [formed by Aβ  only]. Both of these structures can grow to be quite large.   wild-type proteins aggregating into fibrils. Monomers assemble until reaching a critical concentration, for example, a hexamer, as illustrated. Hexamers then aggregate into more complex structures such as protofibrils, and eventually full fibrils; and (c) a model for Aβ  wild-type protein aggregation where monomers join to form dimers, which may grow into protofibrils and eventually fibrils. Protofibrillar intermediates are heterogeneous, metastable aggregates already containing β-sheet regions in the core [38], but retaining some features that are similar to oligomers. The term "protofibril" has varying definitions throughout the literature, where for Aβ the term could refer to structures ranging from 4-11 nm in diameter, up to 200 nm in length, and possibly even longer [39,40]. Protofibrils are considered to be on-pathway during fibrillogenisis, and they could grow larger via monomer addition or merging with other oligomers or protofibrils [41,42]. As protofibrils grow longer, the β-sheet region grows larger. Eventually, a stable and tight β-sheet network is formed in the core by the backbone H-bonds and the hydrophobic interactions of side chains [38]. Figure 2 illustrates the cross-beta structure of fibrils comprised of the prion Sup35. These features are generally associated with protofibrils and fibrils [38]. Another early intermediate thought to play a role in the formation of fibrils is the Aβ protofilament, which can range in diameter from 2.5 nm up to about 6 nm [43], and are 50-100 nm long [44]. Several protofilaments may then merge and form fibrils, which may exhibit a twisted, helical ribbon structure [45] and contain highly-ordered β-sheet regions. A cartoon illustration of an Aβ protofilament is shown in Figure 3. The fibrils may even be composed of several segments with distinct morphologies and varying levels of ordered structure [46]. Additionally, mature fibrils have a diameter ranging from 7-12 nm, and may grower longer than 1 µm [45]. Typically, fibrils are linear, non-branching structures. They contain very large amounts of β-structure, and are generally insoluble. Fibrils can further assemble into bundles [47], and may form plaques outside the neurons. The total assembly process from monomers to fibrils for a simple 1D model is illustrated in Figure 1a, and models for Aβ  and Aβ  fibrils are illustrated in Figure 1b and 1c, respectively. Figure 2. In both (a) and (b), the illustrated cross-β structure is the sequence segment GNNQQNY from the prion Sup35. Carbon atoms are purple or grey/white, oxygen atoms are red, and nitrogen atoms are blue. In (a), cross-β structure is illustrated. Grey arrows represent the back-bone of a β-strand, and the side-chains are shown projecting from the strands. Purple arrows represent the strands residing in the back of the structure. The regions between the strands are referred to as the dry interfaces, whereas just outside of strands are wet interfaces. The fibril axis is indicated by an arrow running through the dry regions between the strands; (b) Side view of the fibril. The H-bonds are formed between red carboxyl groups and blue amide groups from adjacent layers; in (c), a top view of the fibrils shows the interdigitation of two β-sheets, referred to as the steric zipper. Within the steric zipper, water molecules are absent (a red plus sign indicates water). Both images are reprinted from Nelson et al. [48,49].

Statistical Mechanical Approaches to Protein Folding and Aggregation
One of the earlier applications of statistical mechanical methods to protein systems is the Zimm-Bragg model [8][9][10]12] originally developed for the studies of helix-coil transitions in proteins. Although the problems of macromolecular self-assembly that we are dealing with are quite different from conformational changes in proteins, the lesson that can learn from the model is fundamental. The Zimm-Bragg model is like an Ising model, which catches some essential features of a problem, and is solved rigorously using statistical mechanical methods [14]. Extending the ZB model to protein aggregation was first advanced by Oosawa and Kasai [51], Terzi et al. [52], and more recently applied by van Gestel and de Leeuw [25], Schmit et al. [53], and others [54][55][56][57][58][59][60] to the study of macromolecular aggregation. These simple statistical mechanical models may be used to predict the average lengths of protofibrils and fibrils, and the fraction of protein molecules that assume various conformational secondary structures, including sheet, coil, and possibly helix. The main models summarized here focus primarily on computing partition functions for protofibrils and fibrils using simple effective Hamiltonians and transfer matrices.

Partition Function for Helix-Coil Transitions in Proteins
Although protein aggregation is the subject of the article, we use the simplest model, i.e., the Zimm-Bragg model for helix-coil transitions in proteins, to illustrate the power of the statistical mechanical techniques. We first define the coil conformation of a protein residue as the reference state, which means its statistical weight is one [61,62]. The parameter s is the equilibrium constant for a coil residue converting into a helical residue, and relates to the free energy change of adding a helical residue to one end of a helical block, ∆G s = −RT ln s, where R is the gas constant. The nucleation step is to convert a coil residue into a helical residue in a chain of coil residues. The equilibrium constant associated with this event is σs, where ∆G nuc = −RT ln σ is the additional free energy barrier the coil must overcome before converting into the first helical residue that can eventually be part of the helical block. Conventionally, σ and s are referred to as the initiation and propagation parameters. A single helical block in a protein thus has free energy ∆G = ∆G nuc + n h ∆G s , where n h is the number of helical residues in the block.
In direct combinatorial approaches to solving the partition function, Z N for the helix-coil transitions in proteins of length N amino acids is often expressed in terms of the initiation and propagation parameters. Some approximations can be made to simplify the mathematical expression for Z N [62,63]. For example, the simplest model for helix-coil transitions assumes a single helical stretch where all of the residues in the protein are locked into the helix conformation. A more general approach assumes that the conformation of residues could depend on its neighbors and the nucleation and propagation of the chain occurs via a "zipper" mechanism, that is, one single stretch of helix can form along the chain of N residues but may vary in length from one up to N . The partition function for the chain of N residues in the zipper model can be written as: where the term (N − k + 1) is the degeneracy in the number of ways putting k helical residues next to each other along a chain of length N . The Zimm-Bragg (ZB) model for helix-coil transitions in proteins assumes that any number of helical stretches may form along the chain, where residues could be involved in short-ranged interactions with other residues. Each residue could assume either a coil or a helical conformation, thus for nearest-neighbor interactions between residues, there are four possible combinations of a pair of residues. That is, two residues at positions j − 1, and j, respectively, could have states cc, hc, ch, or hh. The jth residue involved in the pair is assigned a weight depending on its conformation, and the conformation of the residue at j − 1. For cc, the weight is one. If two neighboring residues adopt the ch conformations, the helical residue is assigned the weight σs, which corresponds to the nucleation of a helical block. On the other hand, in the original ZB model, the state hc simply has weight one. That is, the growth of a helical block only proceeds in one direction along the chain. Finally, the hh state corresponds with a weight of s at the jth residue. The ZB model is summarized in Table 1. Table 1. Summary of the ZB weights for two residues that are adjacent to each other in a protein.
The partition function for the ZB model, Z N , can be easily computed by using a transfer matrix [64,65]. A transfer matrix is a device that can be used when a system of N units can be decomposed into a subsystem of nearest neighbor, or next nearest neighbor, etc., interactions between all units. Since the residue located at a position j along the chain only depends on the conformation of the residue at chain position j − 1, the transfer matrix factors the partition function for some given energy function E(r) of a system of N units as: where the vectors f | and | i represent the states of the residues at either end of the chain, respectively, λ i and N λ are the eigenvalues, and total number of eigenvalues of T , respectively. Additionally, in the third line of Equation 2, the coefficients c i take into account the effect of the boundary conditions. The notation r = (r 1 , r 2 , ..., r N ) refers to the set of all N spins that represent the states of the residues, where r i denotes the state of the ith residue. The transfer matrix elements represent the probability that a residue occupies a different state from its neighbor. Thus, the transfer matrix used in the ZB model for helix-coil transitions in proteins has the following form: where the column (h, c) is the (j − 1)th state, and the row (h, c) the jth state. As indicated in Equation (2), the matrix T can be readily diagonalized and the partition function written in terms of its eigenvalues. We show the diagonalization as an example in the next section while discussing the helix-coil and sheet-coil transitions in equilibrium protein aggregation. For periodic boundary conditions, the c i 's all equal unity in Equation (2) and the partition function can be written as: where "tr" refers to the trace operation and λ 1 is the largest eigenvalue of Equation (3). Equation (4) is valid only when N is large and becomes exact in the thermodynamic limit when N → ∞. The partition function then reduces to:

Thermodynamic Properties of Proteins
Once the partition function for the protein is known explicitly using the ZB theory, some average quantities can be defined and compared with experiments. The average fraction of residues in a chain of length N that are helical is referred to as the helicity, θ. It can be defined as: where N H is the maximum number of helical residues in the chain. The helicity is akin to measuring the magnetization of spin systems when an external magnetic field is supplied. For the helical protein, the average number of helical segments, v, and the average helical length, L, are also found using Equation (2), Similar averages are calculated for the sheet-coil and helix-sheet-coil systems. In general for a chain of length N , the average number of residues having the property x j is given by: where j could refer to helix, coil, sheet, etc. Thus, the helicity and the average length of helical segments for long chains are now easily computed by inserting Equation (5) into Equations (6) and (8), we are left with: and a similar result can be derived for θ and L when using finite boundary conditions.

Equilibrium Protein Aggregation
In 1961, Oosawa and Kasai constructed a model for equilibrium protein aggregation using ideas from the helix-coil theory that were being developed at the same time. First, in the model the total number of proteins in the system is fixed and is denoted by m tot . Next, the model assumes that a dimer is the smallest aggregate that may form. The chemical reaction for dimer formation represents the nucleation of an aggregate [51,52,66] and can be quantified by an equilibrium constant denoted K eq : where A represents monomer, A k represents the aggregate containing k proteins and is referred to as a k-mer. The nucleation equilibrium constant is often denoted K eq = σs, for example, in Terzi's model [52]. The concentrations for monomers, n 1 , and dimers, n 2 , can also be written as: Once a dimer is formed, then trimer, . . . , k-mer may form by successive addition of a protein to an aggregate. Any of these reactions can be described by the monomer addition mechanism, represented by: where s is the equilibrium constant. Thus, if the equilibrium constant for monomer addition is s, we can write the equilibrium concentration for the k-mer as: Since the total protein mass in the system is conserved, m tot can be written in terms of the concentrations of monomers and aggregates as: Therefore, in the thermodynamic limit N → ∞ the expression can be written as: where the sum converges when sn 1 < 1. If m tot is known, Equation (17) can be solved for the monomer concentration n 1 .

A Generalized Zimm-Bragg Model for Protein Aggregation
A more recent approach to equilibrium peptide assembly introduced by van Gestel and van der Schoot relates the concentrations of protein aggregates to equilibrium partition functions [25].
The partition functions of aggregates are expressed in terms of ZB initiation and propagation parameters and a transfer matrix.
The model describes 1D protein aggregation, where the protein monomer is dominated by coil, sheet, or helical conformations, discussed below, and may participate in short-ranged interactions with other proteins. Hence, the aggregates may exhibit various degrees of conformational order, where helix-coil [24,67], sheet-coil [25,68], or helix-sheet transitions [21,26]

may occur.
This modeling approach is consistent with a recent set of experiments [69], for example, that have shown that oligomers of Aβ(1-40) and Aβ  are dominated by antiparallel β-sheet structures, while their fibrils are mainly characterized by parallel β-sheet structures. Thus major conformational changes may take place somewhere between the oligomer and the fibril formations, and using Ising-like ZB models may be advantageous for studying conformational transitions involved in protein aggregation.
The isolated monomer in the ZB model for aggregation is assumed to be a natively unstructured protein. A "helix" protein is defined if θ helix > θ sheet and θ helix > θ coil , where θ helix can be defined by using Equation (10), or a related equation for sheets. Similar definitions define "sheet" and "coil" proteins. The random coil does not have stable secondary structures. This toy model is suitable for understanding how proteins form fibrils in 1D, however, a word of caution is that, as mentioned above, amyloid formation can be sequence-dependent, for example, Aβ(1-40) and Aβ , have different pathways [27].
All of the aggregate species in a system of volume V are assumed to be in kinetic equilibrium with each other, where interest lies in studying the thermodynamic properties of these aggregates. The system of proteins and aggregates is also assumed to be well mixed, and containing a fixed amount of protein mass given by Equation (16). If the system contains only low concentrations of proteins and aggregates, the solution properties of the aggregates can be calculated by employing a standard ideal gas approximation for the k-mers, where the partition function for the system can be written as: where Z k is the canonical partition function for the k-mer. The number distributions n k per unit volume define the number densities ρ k ≡ n k /V . The relative densities of k-mers can be derived by considering the total free energy density ∆F, which may be written compactly for a system containing N number of proteins, as: which contains an entropy of mixing term as well as the free energy of the aggregate of size N . We can minimize Equation (19) with respect to the total number density ρ T and subject to constraint given in Equation (16), i.e., conservation of mass, which yields for the number densities: where µ is the Lagrange multiplier, and is realized as the chemical potential of a protein, ρ(N ) is just the N th moment of quasi-grand ensemble Ω = N T k=1 ρ(k). As in the 1D model by van Gestel et al. [24,25], the state of an aggregate is directly coupled to the aggregate size distribution.
A generalized ZB model for protein aggregation can now be defined by an effective Hamiltonian. The effective Hamiltonian is used to find a transfer matrix by assuming the interactions between aggregates are described by a nearest-neighbor, Ising-like model, in which the protein could be in any of the two states: a sheet (or helix) or coil conformation. The interactions include the free energy R < 0, which describes the inter-facial tension between adjacent sheet and coil proteins in an aggregate. The parameter P > 0 represents the interaction between two neighboring proteins, where one of the proteins located at position j along the chain is in a sheet conformation. P for sheets is measured relative to the coil interaction energy, which was taken to be zero. Additionally, the free energy K > 0 quantifies a polymerizing interaction between any two monomers along the 1D lattice that does not depend on their respective conformations. The effective Hamiltonian used by van Gestel and others has the following form for N monomers [25,70]: where r = (r 1 , r 2 , ..., r N ) and r i can take on values {1, −1} corresponding to the spin states {↑, ↓} in the Ising model, and to {s, c} in a Zimm-Bragg model for sheet-coil aggregates. With periodic boundary conditions, the partition function for the two-state model can be written as: where e β(N −1)K ≡ k N −1 with β = 1/k B T , and the parameters are redefined as The transfer matrix can be written as: where σ 1 = exp(−2R) and s 1 = exp(P ) are the initiation and propagation parameters for the aggregate system. We note that Equation (23) and Equation (3) yield the same characteristic equation, hence they predict the same thermodynamics results. The difference in the formulation of the folding model versus the aggregation model is that the helical regions of proteins, for example, may only elongate in one direction, while helical aggregates may grow longer at either side. While amyloid fibrils during fibril formation [71,72]. It may also be advantageous to study statistical mechanical models that can describe the interaction of molecules capable of binding to Aβ structure in fibrils, which may inhibit oligomer formation in the early stages of aggregation [73].

Partition Functions for Fibrils
Amyloid formation is generally believed to be dominated by 1D or quasi-1D chains of proteins, which may then bundle into protofibrils and fibrils. It is because of this fact that the transfer matrix formulation in statistical mechanics, if extended successfully, is a powerful technique for the studies of amyloid formation. We focus on this extension in this and the following sections.

Potts Model for 1D Filaments
To include helix, sheet, and coil conformations in a single model, a Potts model [74] for 3-state proteins can be used [26]. The spin variable, r, may now assume values of r = 0 for coil proteins, r = 1 for helical proteins, and r = 2 for sheet proteins. Interactions between proteins are assumed to be with nearest-neighbors only so that a dimensionless Potts model for the aggregate containing N number of proteins can be written as: where the Kronecker delta δ(x, y) equals one if x = y and zero otherwise. Equation (24) is illustrated in Figure 4. Like in Equation (21), the initiation parameters are defined as σ(r j , r j+1 ) ≡ exp(−2R(r j , r j+1 )) and R(r j , r j+1 ) > 0 is the free energy of the interfacial tension between proteins at positions j, j + 1 that are not in the same conformation [24]. Thus, three types of interfaces between neighboring proteins in a generalized model are possible: hc or ch, R(0, 1) = R(1, 0) ≡ R 1 ; sc or cs, R(0, 2) = R(2, 0) ≡ R 2 ; and sh or hs R(1, 2) = R(2, 1) ≡ R 3 . The notation can be simplified by letting σ(1, 0) = σ(0, 1) ≡ σ 1 , σ(0, 2) = σ(2, 0) ≡ σ 2 , and finally The free energies associated with each conformation are listed, as well as the interfacial free energies R 1 , R 2 , and R 3 between helix-coil, sheet-coil, or helix-sheet regions, respectively.
The propagation parameters s 1 and s 2 are associated with the free energies P 1 and P 2 that refer to the interaction between the ith protein that is helix or sheet, respectively, and the nearest neighbor protein at location i + 1. The coil protein interaction energy is assumed to be zero so that it may serve as a reference state. The transfer matrix can then be written as and the partition function for N > 2 proteins in the aggregate can be calculated by diagonalizing Equation (25) and plugging into Equation (2). The Ising-like ZB model for sheet-coil (or helix-coil) transitions in aggregates, Equation (21), can be recovered by writing the q = 2 version of the effective Hamiltonian given in Equation (24). By choosing internal states such that r i = −1 is the coil state and r i = +1 is the helix state, and by making the substitution δ(r i , r j ) = 1 2 (1 + r i r j ), the effective Hamiltonian and corresponding transfer matrix for 1D, two-state helix-coil is recovered [24]. Next, the Hamiltonian for the 1D model given by Equation (24) can be applied to describe the interactions between proteins in aggregates on quasi-1D lattices.

Simple Model for Fibrils
To describe the fibrils, which may contain several filaments that can be described by Equation (24), various models could be used. For example, two identical filaments could align in register, and all of the proteins could be in the sheet conformation already, or all coil, or some combination of sheet, coil, and even helix proteins. A simple Hamiltonian for the all-sheet case can be written as [25,55,68]: where the free energy B > 0 describes a lateral binding interaction between two sheet proteins on different filaments, L y refers to the the number of filaments, and L x is the length of each filament. Additionally, in our approach K was the polymerizing interaction between two adjacent proteins in a fibril, and P 1 was free energy of an interaction between a sheet protein and one of its neighbors.
Other effective Hamiltonians could also be written to describe the case where the filaments are not aligned in register with each other, as well as cases where the protein conformation may also play a role in the assembly of the filaments into full fibrils [55,68].

Quasi-1D Models for Aggregates
In addition to conformational changes in filaments and fibrils, another complication is that the kinetic pathways to the formation of fibrils seem to be sequence-dependent [27]. The Aβ(1-40) isoform solution is abundant in dimers, then trimers, tetramers, . . . , in decreasing order. However, the Aβ(1-42) isoform is more abundant in hexamers and pentamers than in dimers and trimers [4,27]. These facts seem to be consistent with recent experiments [38,69,75], which indicate that the Aβ(1-40) dimer is particularly stable and contributes to protofibril formation. On the other hand, circular hexamers seem to play a role in the protofibril formation of Aβ   [4,27].
We can model the equilibrium aggregates of Aβ(1-40) and Aβ(1-42) proteins by using finite strips of a two-dimensional N × N square lattice. In Figure 5 and Figure 6a, two identical 1D lattices stacked in-register, that is, a strip lattice of width two, are used to represent an Aβ(1-40) aggregate. For Aβ(1-40), we assume that the smallest equilibrium aggregate is the critical nucleus, which in this case is taken to be the dimer, while the smallest aggregate for Aβ  is the hexamer, as illustrated in Figure 7 for a strip lattice of width six. The position of a vertex within the strip is specified by coordinates (i, j), where i is the position along the x-axis of length L x vertices and j is the position along the y-axis of width L y vertices. The total number of vertices is N T = L x L y . In Figure 5a, strips of spin variables s j i in the y-axis are referred to as L y -mer's, where s j i = −1 or +1 for Ising-type models and 0, 1, or 2 for Potts models. The critical nucleus can be represented by a column of L y proteins on the strip lattice. The proteins in the nucleus could also participate in inter-protein interactions other than the polymerizing interaction K in the y-direction. These interactions can be described by using the sheet and helix interactions from the filament model, plus the free energy introduced above, B > 0, that quantifies the inter-filament interactions between two sheet proteins. The interactions between the proteins in these aggregates is modeled similarly to the 2-helix chain model for proteins proposed by Skolnick [76] and others [59,63,77,78], which use Zimm-Bragg or Lifson-Roig (LR) parameters to quantify the inter-chain interactions between residues in independent chains. When the inter-chain interactions between two helical residues are made zero, for example, the partition function reduces to a direct product of Zimm-Bragg [76] (or LR [63,77]) transfer matrices. Since the model for aggregates proposed here uses a strip of a finite 2D lattice, as illustrated in Figure 6a, the two-protein case studied by Skolnick can be considered as a special case when considering protein folding instead of aggregation. The lessons learned from the folding models can guide us in constructing aggregation models. Figure 7. (a) A strip lattice with periodic boundary conditions in the y-axis is used to model the nucleus for Aβ(1-42) as a hexamer. The strip has L y monomers per column in general. B (blue, dotted lines) is the free energy contribution from the interaction between proteins s j i , s j+1 i that are both sheet along the y-axis; (b) Oligomers aggregate along the x-axis with a total of L x sites, where L x → ∞ is the thermodynamic limit.
As mentioned, the nucleus is the smallest equilibrium aggregate in our formulation. The next smallest aggregate occupies the first two columns of the strip lattice, and contain 2L y number of proteins, then 4L y number of proteins, and so on. The interactions between L y -mers along the x-axis can be described by generalizing the effective Hamiltonian for the 1D aggregation model, that is, P 1 > 0 and P 2 > 0 represent the interaction between helical and sheet proteins, respectively, and their nearest-neighbors in the aggregate. The total effective Hamiltonian for aggregates on a strip lattice with boundary conditions can be written as: where H f il is given by Equation (24) upon substituting r i → s j i and where s j i can assume the values of 0, 1, 2 for coil, helix, or sheet monomers, respectively. Additionally, B is the lateral binding interaction between two proteins from different filaments. The third term involving K are the polymerizing interactions between proteins in both the x and y directions on the quasi-1D lattice. We assumed any polymerizing interactions between proteins on a quasi-1D lattice were equal in magnitude in order to keep the number of parameters used in the model to a minumum. It is similar to the parameter K in Equation (21), which took into account the polymerizing interaction between two adjacent proteins on the 1D lattice. For the case q = 2 and L y = 2 of Equation (27), the Hamiltonian for sheet-coil protofibrils can be explicitly written as: where 2D refers to aggregates on a strip lattice. The transfer matrix can then be written as: where b ≡ exp(B). As mentioned, in the limit when the inter-filament interactions between two sheet proteins B → 0, the transfer matrix given by Equation (29) decomposes into a direct product of transfer matrices given by Equation (23) for 1D filaments of length N along the x-axis [25]. Moreover, the limit in which the sheet-coil interfacial interactions R → 0 (or σ 1 → 1) yields independent strips parallel to the y-axis. The 2 Ly × 2 Ly transfer matrix is also symmetric with respect to the sheet-coil interfacial interaction R. This fact is analogous to the 1D case, where the transfer matrix was symmetric in σ 1 , σ 2 and σ 3 . The partition function for aggregates on the strip lattice can be calculated by plugging the eigenvalues, λ 2D,i , of Equation (29) into Equation (2) and specifying boundary conditions. The result is: where, again, c i are determined by boundary conditions and k was defined in Equation (22). Using the ZB formalism, we can also model fibrils in the x-, yand z-directions using a quasi-1D lattice in 3D. For example, two or more filaments could join to form a protofibril or fibril, as illustrated in Figure 6 for two simple geometric configurations. For simplicity, we study the case where two 2D aggregates represented by strip lattices, as depicted in Figure 6b, are the same geometrical shape and are stacked one on top the other, in-register. The spin variable associated with one of the aggregates is denoted by s, while the spin variables associated to the second aggregate is denoted by t. The effective Hamiltonian for the fibril model in Figure 6b, referred to as the "cube" model, may be written using the strip model for aggregates as: where we assume that the interaction between any two sheet proteins from adjacent aggregates is also quantified by the free energy B > 0. Additionally, the aggregates can now polymerize in any direction.
To help keep the number of parameters used in the model to a minimum, we assume the polymerizing interactions between two adjacent proteins in the z-direction has the same strength as the polymerizing interactions, K, in the x and y directions. The corresponding transfer matrices for each model for fibrils are found just as they were for the 1D and strip models discussed earlier. The result is: In general, the transfer matrix for the 2D model has dimension q Ly × q Ly , while the 3D model has dimension q 2Ly × q 2Ly , where we assumed that the two protofibrils composing the fibril contain L y number of filaments. The transfer matrices for both the 2D and 3D models are illustrated in Figure 6a and Figure 6b, respectively. Since the Potts model is used, q is 2 for sheet-coil, helix-coil, or helix-sheet models, and 3 for the helix-sheet-coil model.

Dilute Thermodynamic Averages
To compare with experiments, some average properties of the dilute system of monomers and aggregates can be defined. The total number density for the strip or cube model for fibrils can be written as: where Z can be given by Equation (30) or (32) for 2D or 3D aggregates, respectively, and L is the total number of proteins in aggregates. In the 2D case, L = L x L y . The average fraction that an aggregate is helix (i = 1) or sheet (i = 2) can be defined as: where the x-axis is the axis of propagation of the aggregate, and θ i can be calculated for any of the aggregate species discussed earlier by computing Equation (9). Equation (34) can be used as a fit function for CD spectral data points. The average degree of polymerization of an aggregate growing in the x-direction can also be defined as: and is directly related to the length of the fibrils. The expressions for θ 1 and L can be obtained for systems of α-synuclein (αS) and Aβ  aggregates. The α-synuclein fibril is modeled by placing the proteins in the aggregates onto the L y=4 strip lattice, as illustrated in Figure 6a. Thus, the average length of these fibrils is then: which can be used to fit the AFM measurements of the average lengths of the fibrils. This relation also holds for the average length of a fibril described by the cube lattice model as depicted in Figure 6b for the Aβ(1-40) fibrils. As an example, we calculate these quantities explicitly for a system of equilibrium Aβ(1-40) aggregates using a sheet-coil model. Specifically, for a system where 1D filaments, L y = 2 strip aggregates, or 3D cube aggregates could be present at equilibrium, we first define: where the fugacity z ≡ exp(βµ) and in each sum λ k is the kth eigenvalue of the 1D, 2D, or 3D transfer matrix for filament, strip, or cube models, respectively. The Aβ  aggregates considered here at equilibrium are illustrated in Figure 8. Additionally, x k is the kth term of the expression i|f , where |i and |f are the specified boundary conditions. The sums computed converge only if kzλ j < 1 for all j. Details on boundary conditions for ZB-type models can be found in [67]. By using Equation (16), φ = m tot /V can also be written explicitly for Aβ  aggregates as: Now the average lengths of the aggregates can be computed. Next, the average fraction of the aggregates that are sheet, θ 2 , is calculated for the Aβ(1-40) model. By plugging Equation (6) for filament, strip, and cube aggregates into Equation (34), θ 2 can be written as: where s 2 is the sheet propagation parameter. The procedure for finding θ 2 is quite general, and works for all of the transfer matrices that we have considered in this model. Figure 8. Partial lists of chemical species that may exist in dynamic equilibrium with fibrils for Aβ(1-40) (top) and α-synuclein (bottom). In the Aβ model, the different types of aggregates that could be present at equilibrium are 1D filaments, strips of length L y = 2 that represent protofibrils, and 3D cubes that represent fibrils. The cubes are composed of two identical proto-fibrils stacked in-register. In the model for α-synuclein aggregates, L y = 1, 2, 3, and 4 strip lattices are used to describe the aggregates at equilibrium, with the L y = 4 strip lattice representing the fibril. For both Aβ(1-40) and α-synuclein, we assumed n c = 2.

Comparison to Experiment
In this subsection, our ZB-like model predictions are compared to the experimental results for the CD spectra of Aβ(1-40) fibrils [52] and the AFM measurements of the lengths of α-synuclein fibrils [68]. The CD and the AFM measurements were made at various initial mass concentrations of each protein, when the fibrils had reached a steady state. For the fit of the CD data, we used as our fit function the total fractional amount of sheet structure in all of the aggregates, θ 2 Aβ , given by Equation (39). The proteins at the boundaries of aggregates could be coil or sheet for 1D, 2D, and 3D lattices. The fit and the values for P 2 , the sheet interaction free energy, K, the free energy describing the polymerization of the aggregate in any direction, R 2 , the sheet-coil interfacial free energy, and B, and lateral binding free energy for the q = 2 sheet-coil model are given in Figure 9a. The fit parameters are then used to predict the average length of the fibrils in the system, L , which illustrated in Figure 9b. For the α-synuclein model, where L y = 1, 2, 3 and 4 strip aggregates could be present, we fit the L y = 4 contribution from Equation (36) to the AFM average length data for the fibrils. The α-synuclein aggregates at equilibrium are illustrated in Figure 8. The fit is illustrated in Figure 9c. The fit parameters are then used to predict the total fraction of aggregates that are sheet, which can be written for the α-synuclein model as: where each contribution in Equation (40) can be calculated using transfer matrices derived from Equation (27) for L y = 1, ..., 4. Equation (40) is illustrated in Figure 9d. The model predicts the average length data pretty well (Notice the log scales used in the plot. At the low coverage, the predicted points fall actually within the error bars), as we should have expected because the other variations of the Ising-ZB model fit the data well [25,68]. The resulting predictions for θ 2 αS illustrate that the concentration at which the fibril concentration takes off is around 15 µM, again, as we should have expected [68]. The fit of θ 2 Aβ to the CD data predicts that the fibrils are held together tightly due to the relatively large value of the sheet-interaction between proteins in the aggregates, P 1 , and to a lesser extent on the binding between filaments as quantified by B. The fitted value for the sheet-coil interface free energy, R 2 , indicates that the interfacial tension between sheet and coil regions in the aggregates is modest. However, the fibril concentrations do not really increase from zero until nearly 100 µM according to the model predictions, but β-rich fibrils have been observed at lower concentrations as seen in Figure 9a. The fit could be improved by considering other types of models for the Aβ fibrils and as well as different boundary conditions. The model predictions for the strip model of α-synuclein fibrils seem to agree with the experimentally determined average lengths as illustrated in Figure 9c where the boundary conditions were set so that the ends of fibrils could be sheet or coil proteins. Additionally, as illustrated in Figure 9d, the model for synuclein fibrils predicts that the sheet-coil transition of proteins in fibrils largely drives the polymerization process, where ρ f ib /φ and θ 2 Ly=4 give nearly the same result for the concentrations used in the AFM experiments.
The fits of the AFM data done by van Raaij et al. [68] and Schmit et al. [53] needed only 2 parameters, compared to 3 in the present model. When compared with van Raaij's fit of the AFM data, our model predicts that the probability that fibrils contain sheet structure is high once overcoming a sheet-coil free energy barrier R 2 , whereas in van Raaij's model prediction, the free energy barrier between adjacent sheet and coil proteins in the aggregates does not seem to be present. A finite contribution from R 2 means the fibrils will have longer stretches of sheet content when compared to cases when R 2 is closer to zero when there is little or no penalty between coil to sheet regions in the fibrils. Additionally, our model predicts that the inter-filament interactions, B, are slightly weaker than the interactions between sheets along the axis of growth of the fibril, whereas van Raaij's fit is the other way around: the inter-filament interactions are stronger than those between proteins along the fibril axis of growth.
When fitting the Aβ(1-40) CD data, our model predicts that the value of the polymerizing interaction between proteins on a quasi-ID lattice, K, is small and could be due to modeling uncertainty, thus the number of parameters needed to fit the CD data could be less. The fact that not all adjustable parameters were required to fit CD and AFM data suggests that the model Hamiltonians introduced throughout the paper may be simplified to describe only the relevant interactions quantified by the non-zero fit parameters. It may also imply that the fibrils are mainly held together by a few types of energetic interactions, for example, the inter-filament interaction B for Aβ(1-40) had a finite contribution in the Hamiltonian. Other interactions described by the Hamiltonians could be non-existent or very small. More detailed experimental results are needed to discern the correctness of the models in predicting the values of interaction energy parameters, which in turn describe the dominant interactions within the fibrils.
As mentioned, the data sets were fit using various boundary conditions, and the open boundary case (proteins could adopt either conformation at the boundaries) was found to be the best choice for fitting for the average lengths of the α-synuclein fibrils. The CD fits could not be shown to be dependent on certain boundary conditions since there is currently no AFM measurements of the fibrils to compliment the CD data. This means we could fit the CD data using most choices for boundary conditions, including the case where all proteins at the ends of fibrils are in the coil conformation. However, for some choices of boundary conditions [25,67], the corresponding average length predictions yielded unreasonable lengths (not shown) for protein aggregates in vivo.

Grand Canonical Approach
The models for fibrils discussed so far do not take into account interactions between protein and solvent, or some free energy that would be associated with nucleus formation. The ZB model for aggregation can be extended to take into account these phenomena by using a grand-canonical model. We summarize several main differences between canonical and grand-canonical approaches: (1) in the grand canonical model, aggregates of all sizes are included; (2) an aggregate phase and solution phase are in equilibrium; (3) chemical potentials can be used relating the solution phase as well as the aggregate phase [26].
The solution phase is defined by specifying the chemical potential for protein monomers in the solution can be written [79][80][81] as: where the subscript "S" stands for solution, µ ST and µ SR are the free energy contributions arising from the translational and rotational motion of monomers moving in solution, respectively, and c is the concentration of monomers in solution. The aggregate phase is defined by specifying the chemical potential of the aggregates, µ agg , by assuming a crystalline approximation so that µ agg can be written as [82]: where "P" stands for polymers of proteins. µ P C is the free energy contribution arising from the contact interactions between proteins in aggregates, which may vary for different monomer organizations in the aggregates. We assume this term also includes the conformational entropy of the backbone and side-chains. The term µ P V is the free energy arising from the proteins vibrating about their equilibrium positions, but not molecular internal vibrations within the proteins [79,82]. When the phases are at equilibrium, the chemical potentials for each phase are equal: With the simple statistical mechanical model summarized in the sections to follow, we can relate the chemical potential contribution from the protein interactions in aggregates, µ P C , to the experimental concentration of protein in solution via Equation (43).
As a first step in generalizing the canonical effective-Hamiltonian models presented earlier, the free energy A is introduced to quantify the entropic penalty needed to nucleate the aggregate, i.e., the first column of a strip lattice that contains protein aggregates. This free energy may also be viewed as a boundary between proteins and solvent. The aggregate phase then assumes strip or cube lattices may be occupied by aggregates and any other species, including solvent clusters.
To write down an effective Hamiltonian that can include the free energy A, the 1D or quasi-1D lattices used in constructing fibrils can be generalized to allow solvent clusters to occupy the lattice sites. For example, in Figure 10, a square represents a solvent cluster, whereas circle represents protein. Both solvent and proteins can occupy sites along 1D or quasi-1D lattices. By introducing a lattice gas model into the aggregate phase, a Potts Hamiltonian for the 1D lattice that quantifies the interactions between helix, sheet, or coil proteins and solvent can be written as [26]: −βH nc ps = − where the lattice-gas variable n i = 1 refers to a protein occupied lattice site, and n i = 0 a solvent occupied site. Additionally, "pp" in −βH pp refers to "protein-protein" interactions and "ps" in −βH nc ps refers to "protein-solvent" interactions. The term χ(n i , n i+nc ) = 1 − δ(n i , n i+nc ) ensuring that there is solvent at site i and a protein at i + n c , or vice-versa. Figure 10. Summary of protein conformation energies. A site could be occupied with a solvent cluster, denoted by n = 0 (square), or a protein, n = 1 (circles). Proteins may assume a particular conformation (sheet, black/solid circle; coil, white circle). A dilute q = 2 Potts model for sheet-coil conformations is shown, where n c = 1 and the free energies P 1 , K, R, and A are illustrated.
Since the number of proteins on the lattice can fluctuate, this description of protein aggregation is described by using the grand canonical ensemble. The lattice-gas formalism, i.e., Equation (44), is able to describe a variety of elongation mechanisms including merging and fracturing of aggregates of different sizes along the 1D lattice. The partition function can be written as: where βµ P C is the dimensionless chemical potential arising from the contact and interfacial interactions between proteins in aggregates, and where the sum is performed over both spin and lattice-gas variables. Just like in the canonical models, Q may be solved for exactly by a transfer matrix T . A simple example illustrating T for the case n c = 1 in a sheet-coil (t i = −1 for coil, t i = 1 for sheet) system can be written as: where s 1 , σ 1 , and k were defined earlier and α ≡ exp(−2A) is a new Zimm-Bragg-like parameter. Additionally, the fugacity is now defined as z ≡ exp(βµ P C ).
The inter-filament interactions between two 1D filaments are treated using the same methodology introduced in earlier sections. In general, the Hamiltonian for an L x × L y strip lattice that includes inter-filament interactions can be written using the 1D Hamiltonian, Equation (44), by changing the spin and lattice-gas variables t i → t j i and n i → n j i , respectively, as: where H f il (j) refers to the jth filament. For Aβ(1-40) L y = 2, as illustrated in Figure 5b,c. The parameter F quantifies the interaction energy between two sheet-linked proteins from adjacent filaments, and plays the same role as the free energy B in earlier models for fibrils in this article. In our treatment F > 0, the proto-fibrils and fibrils are more stable than single filaments. Since nucleation cannot in reality occur in 1D, we consider a similar model for aggregates that positions the nucleus along the y-axis, as shown in Figure 6a and Figure 11. From this point of view the orientations of proteins in the nucleus are perpendicular to the direction of propagation (x-axis) of the fibrils, and the nucleus is now a multi-layer, 1D aggregate. The nuclei may assemble into proto-fibrils that grow longer on the quasi-1D lattice. An effective Hamiltonian for protein aggregation, including the quasi-1D nucleus, can be written: where −βH pp (j) is given by Equation (45) after substituting t i → t j i and n i → n j i . In the y-direction we write analogous interactions, −βH y , similar to those in the x-direction. Also included in the y-direction is the nucleus term containing the parameter A, which has the same meaning of surface energy as before. The effective Hamiltonians given by Equations (49) and (50) are the most general forms of fibrils that we have considered so far. Figure 11. Proteins or solvent clusters may occupy lattice sites, where the front-view (y-z plane) of an aggregate of Aβ  proteins is shown along with the interactions between proteins and solvent clusters. The n c = 2 nucleus is represented by dashed-dotted lines (free energy A denoting the nucleation). Dotted and solid lines illustrate interactions between sheet proteins. Double solid lines illustrate a protein-solvent interface. Dashed (blue) lines have no meaning.
For either description of fibrils (model A or B), the total number of proteins on a strip lattice is then Ly j=1 n j i . The grand partition function can be written as: where the sums over {t}, {n} are for all i and j, and A, B refers to the effective Hamiltonians given by Equation (49) or (50), respectively. For periodic boundary conditions, Equation (53) can be solved as strip is the partition function for the lattice-gas model (A or B). Just as in Subsection 3.1, in the thermodynamic limit N T → ∞: is the largest eigenvalue of T A(B) strip . In general, the dimension of the transfer matrix T A strip is (q + 1) ncLy × (q + 1) ncLy and has (q + 1) ncLy number of eigenvalues, whereas the transfer matrix T B strip is (q + 1) Ly × (q + 1) Ly and has (q + 1) Ly number of eigenvalues.
To compare with experiments, we can define quantities similar to Equations (34) and (36), and others. For example, in the grand canonical ensemble, the average number of proteins on the lattice, N p , referred to as the occupation of the lattice, the number of proteins in filaments, ψ , the average number of filaments, γ , and the average number of sheet segments, θ , can be written as: respectively. Other quantities may also be defined including the average lengths of filaments, and the average length of sheet stretches in aggregates [26,83].

Comparison to Experiment
We solve for µ P C for Aβ(1-40) from Equation (43) in terms of µ ST , µ SR , µ P V , and experimental concentration c. Then, µ P C is plugged into Equation (53), and Equation (57) and other thermodynamic quantities can be calculated. For Aβ(1-40), we have µ ST + µ SR ≈ −29 kcal/mol [79,81]. In reference [79], µ P V for hemoglobin was found to be approximately 3 4 (µ ST +µ SR ). Mutated hemoglobin does polymerize as amyloid, but amyloid proteins usually may be natively unstructured, unlike hemoglobin. We nevertheless use a similar result for µ P V for Aβ(1-40) and Curli. Equation (57) divided by Equation (55), θ / N p , the β-sheet fraction, is used as our fitting function. The results are plotted in Figure 12a for Aβ(1-40) fibrils, and Figure 12b for the Curli fibrils. The fit yields reasonable free energies at room temperature for the Aβ(1-40) fibrils, P 1 ≈ K ≈ A ≈ 0 kcal/mol, R 1 = 0.35 kcal/mol, and F = 16.4 kcal/mol. For the Curli fibrils, we found P 1 = 7.26 kcal/mol, K = 2.2 kcal/mol, R 1 ≈ 0 kcal/mol, and A = 1.2 kcal/mol [26]. Clearly for the experiment involving Aβ(1-40) fibrils, the grand canonical approach to modeling fibrils does a better job than earlier canonical approaches. The grand-canonical model also suggests that the Aβ fibrils are more strongly held together by inter-filament interactions when compared to the fit from the canonical model, and also that the penalty in going from sheet to coil regions in the aggregates is very small and could be due to modeling uncertainty. The minimum number of parameters needed to fit the CD data is also the same number when compared with the van Raaij and Schmidt's models, and if the value of R 1 in the Aβ fit is taken to be within modeling uncertainty, then only one parameter is needed to fit the CD data.  [52] involving Aβ  aggregates; (b) The fraction of sheet proteins in Curli fibrils is fitted to the scaled results of the Hammer et al. experiment [84]. In (a), the fit parameters were P 1 ≈ K ≈ A ≈ 0 kcal/mol, R 1 = 0.35 kcal/mol, and F = 16.4 kcal/mol.; while in (b) we have P 1 = 7.26 kcal/mol, K = 2.2 kcal/mol, R 1 ≈ 0 kcal/mol, and A = 1.2 kcal/mol [26]. In (a) we used case B of the strip models with n c = 2 whereas in (b) we used the 1D model with n c = 2 for aggregation. In both cases q = 2.

Conclusions
By focusing on the aggregation of proteins in forming oligomers, protofibrils, and fibrils, and their relations to neurodegenerative diseases, statistical mechanical approaches to protein aggregation have been developed. We have made a general summary of the field, presenting recent formulations of the ZB model based on the canonical, as well as the grand canonical, approaches to the amyloid formation processes. Some results are presented to show that these models can be used to interpret experimental observations as well as to provide phase diagrams [26] showing the parameter dependence of the β-sheet dominating regions.
More experimental data like the CD results of Terzi et al. [52] and the AFM measurements performed by van Raaij et al. [68] would help validate the ZB approach. For example, the ZB model for protein folding has been used to classify all the amino acids based on their propensity to fold from coil to helix [12]. Similar classification schemes could potentially be devised for the many proteins that can aggregate to form fibrils if a much larger collection of experimental data (like the CD and AFM results) were available.
Of course, a statistical mechanical approach to protein aggregation has serious limitations. It can only be used to study the equilibrium properties of the systems, not the rates of the processes involved nor the transient behaviors, such as quasi-equilibrium or kinetic trapping. Furthermore, the available experimental data that we can compare our theories with are so far extremely limited. Therefore, statistical mechanical models are not a tool for predicting assembly pathways. However, statistical mechanics can be used to show that experimental observations are consistent with the predictions of a certain route of aggregation, that is, given a route of aggregation and the associated effective free energy, statistical mechanical models can predict equilibrium distributions of oligomers, protofibrils, fibrils, etc., which can then be compared to observed data.
The pathway prediction function is better achieved by using kinetic models [7,41,42,85,86], or better still, by molecular dynamics simulations [1,2,32,33,87,88]. Unfortunately, the latter are highly restricted by the system size and length of time that molecular dynamics can be used to simulate.
On the other hand, protein aggregation, as discussed earlier, is a complex process spanning many levels of structure and many chemical species. Thus, at present, the kinetic models or coarse-graining models may be better tools for the purpose of predicting pathways. Moreover, many more experimental data are available in the literature for non-equilibrium or kinetic studies. It would be interesting, for example, to develop a kinetic approach based on statistical mechanics, similar to a kinetic Ising model derived from an Ising model. The kinetic study of protein aggregation is a rich field of investigation and of great current interest. We hope to be able to report progress in the non-equilibrium studies of protein aggregation in future work.