This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.105.074906v1 90/7/2510    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Park, J. Articles by Hwang, W. Search for Related Content PubMed PubMed Citation Articles by Park, J. Articles by Hwang, W.

Atomistic Simulation Approach to a Continuum Description of Self-Assembled ß-Sheet Filaments

Jiyong Park ^*, Byungnam Kahng ^*, Roger D. Kamm   and Wonmuk Hwang  

^* School of Physics and Center for Theoretical Physics, Seoul National University, Seoul, Korea; Biological Engineering Division, Center for Biomedical Engineering, and Department of Mechanical Engineering, and Biological Engineering Division and Center for Biomedical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts; and Department of Biomedical Engineering, Texas A&M University, College Station, Texas

Correspondence: Address reprint requests to W. Hwang, Tel.: 979-458-0178; E-mail: hwm{at}tamu.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX: BENDING STIFFNESS OF... ACKNOWLEDGEMENTS REFERENCES

 
We investigated the supramolecular structure and continuum mechanical properties of a ß-sheet nanofiber comprised of a self-assembling peptide ac-[RARADADA]₂-am using computer simulations. The supramolecular structure was determined by constructing candidate filaments with dimensions compatible with those observed in atomic force microscopy and selecting the most stable ones after running molecular dynamics simulations on each of them. Four structures with different backbone hydrogen-bonding patterns were identified to be similarly stable. We then quantified the continuum mechanical properties of these identified structures by running three independent simulations: thermal motion analysis, normal mode analysis, and steered molecular dynamics. Within the range of deformations investigated, the filament showed linear elasticity in transverse directions with an estimated persistence length of 1.2–4.8 µm. Although side-chain interactions govern the propensity and energetics of filament self-assembly, we found that backbone hydrogen-bonding interactions are the primary determinant of filament elasticity, as demonstrated by its effective thickness, which is smaller than that estimated by atomic force microscopy or from the molecular geometry, as well as by the similar bending stiffness of a model filament without charged side chains. The generality of our approach suggests that it should be applicable to developing continuum elastic ribbon models of other ß-sheet filaments and amyloid fibrils.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX: BENDING STIFFNESS OF... ACKNOWLEDGEMENTS REFERENCES

 
Peptide self-assembly has been recognized as a new fabrication modality to generate biomaterials with defined physicochemical properties (1–5). In particular, a class of ß-sheet forming peptides have been designed that can construct filamentous structures (6–9). They are typically 8–16 amino acids in length, and have an alternating sequence of hydrophobic and hydrophilic residues that form a bilayer of ß-sheet tapes with hydrophobic side chains between the tapes. The hydrogel made up of these filaments shows promise as a three-dimensional cell culture matrix or as a tissue engineering scaffold. Due to the short sequence of the constituent peptides, they are easy to synthesize and manipulate to achieve the desired functionality of the hydrogel; the molecular composition and stiffness of the hydrogel can be controlled by varying peptide sequence and its concentration (10), and it is also possible to decorate the peptide with various ligands for further functionalization (11). The ß-sheet peptide hydrogels have been used as scaffolds for neurons (12), neural progenitor cells (13), chondrocytes (14), and endothelial cells (15).

Another important aspect of ß-sheet peptide self-assembly is its similarity to amyloid fibrils found in various protein misfolding diseases (16–19). Regardless of the protein sequence or length, most amyloid fibrils have similar diameters of 10–20 nm (20) and share the cross-ß structure, where the ß-sheet runs perpendicular to the fibril axis (21), essentially the same as those formed by self-assembling peptides. Although the early oligomers in the assembly process rather than the fibrils seem to be the toxic species in neurodegenerative diseases (22,23), in other classes of systemic amyloidoses, the accumulation of fibrils by itself can be symptomatic through compression of blood vessels and adjacent tissues (24,25). In addition, aggregation of amyloid precursor protein (26) or tau (27) can impede axonal transport (28). Peptide self-assembly is thus a good model system for amyloid fibril formation.

In developing the self-assembling peptide hydrogel as a biomaterial, it is useful to have some degree of control over its mechanical as well as chemical properties. Recent evidence suggests the importance of the mechanical environment in determining cell behavior. Substrate dimension was found to affect the polarization of cells (29) or the nature of cell-substrate contacts (30). Moreover, rigidity of the substrate affects the shape and migratory behavior of cells (31–33). Since hydrogels are comprised of a network of peptide filaments, it is useful, as a first step, to study the mechanical properties of a single filament. Although few attempts have been made to directly measure the elastic properties of a single peptide filament, Leon et al. (34) used an elastic strut model to deduce single filament properties from macroscopic rheometry data. Since such an approach is model-dependent, it is still desirable to develop methods to directly probe single filament properties.

In the case of cytoskeletal filaments, various experimental and computational approaches have been used to quantify their elastic properties. Bending stiffness has been estimated by monitoring the thermal motion of fluorescently labeled beads on F-actin (35) or a microtubule clamped at one end (36). A similar method was employed to measure torsional rigidity of F-actin (37). There also have been computational approaches, mainly based on normal mode analysis with simplified force fields, to calculate the elastic properties (38) or to monitor long wavelength motions (39) of F-actin.

Unbranched polymers inherently exhibit a large length/thickness ratio, hence can be described globally as linear chains characterized structurally by a persistence length (40). In the particular case of biologically active filaments, however, molecular details are also important since interactions between different filaments or between filaments and other proteins (e.g., receptors or cross-linking proteins) are often local and highly specific. A model of the biofilament that captures both the atomistic and continuum properties would thus be very useful.

In this study, atomistic simulations are used to investigate the structure and elastic properties of filaments comprised of the self-assembling peptide RAD16II, with the sequence

with N- and C-termini acetylated and amidated (Fig. 1 a). The selection of this particular sequence was initially motivated by its similarity to RGD, an integrin-binding epitope in the extracellular matrix (41). Subsequently, hydrogels formed from RAD16II have been used as a substrate for various cell culture studies (12,14,15). Here, we first determine the possible supramolecular structures of the filament, then characterize its elastic properties using three different simulation methods: thermal motion analysis (TMA), normal mode analysis (NMA), and steered molecular dynamics (SMD).

View larger version (86K): [in this window] [in a new window]   FIGURE 1  (a) Molecular structure of RAD16II. Top part of the molecule (ARG and ASP) is hydrophilic, while the bottom part (ALA) is hydrophobic. We used VMD (65,66) for all molecular visualizations. (b) AFM image of the network of the RAD16II filaments. The scan size is 1 µm. (Courtesy of W. Jeong.)

Here, we first analyze the supramolecular structure of the RAD16II fiber following the approach of Hwang et al. (48); out of all possible combinations of ß-sheet filaments, we select the most stable ones by running molecular dynamics (MD) simulations. Within numerical accuracy, four filament structures that differ only in backbone hydrogen-bonding register are found to be similarly stable. Each configuration exhibits linear elastic behavior within the limits tested by the simulations, and possess Young's moduli in the range 5.5–9.7 GPa for deformations in the transverse directions. The filaments are anisotropic, however, in the sense that they are nearly inextensible in the axial direction. The effective width (5.7 nm) and thickness (1.4 nm) of the corresponding continuum ribbon has an aspect ratio of 4, larger than that estimated from atomic force microscopy (AFM) experiments or from geometrical considerations of the ß-sheet bilayer (6-nm wide and 2-nm thick). Furthermore, simulations of the filament without charged side chains yielded bending stiffness similar to that of RAD16II. These are attributed to the backbone hydrogen-bonding interactions that dominate the elasticity of the filament.

Our findings suggest that although side-chain interactions are important in determining the registry of peptides and equilibrium stability of a ß-sheet filament, they have little influence on filament elasticity, so that use of a generic continuum ribbon model is justified in further developing a network model comprised of individual filaments.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX: BENDING STIFFNESS OF... ACKNOWLEDGEMENTS REFERENCES

 
Simulation methods
For simulations, we used CHARMM (49) version 29 with the PARAM19 force field. For the MD simulations used to determine filament structure, and for TMA and SMD, the solvation effect was incorporated by using the analytic continuum electrostatics (ACE2) module in CHARMM (50–52). For NMA, the distance-dependent dielectric constant method (RDIE) was used instead of ACE2. The choice of different filament lengths below was mainly determined by the computational loads that vary among different simulations.

Identification of the most stable structures
We constructed the candidate filament structures containing 60 peptides each, as detailed in the next section. Each filament was initially energy-minimized with 200 steps of the steepest-descent method, followed by 6000 steps of the adapted-basis Newton-Raphson (ABNR) method. The system was then heated from 0 K to 300 K in 60 ps, and equilibrated for 300 ps by rescaling the velocities to keep the temperature in the range 300 ± 10 K. The final production run without velocity rescaling was performed for 100 ps and the coordinate trajectory was saved every 0.8 ps. The long equilibration period was necessary to ensure relaxation of the initial structure that was constructed by aligning peptides in extended conformations. In both MD simulations (TMA and SMD), the lengths of the bonds connecting hydrogens and heavy atoms were fixed by using the SHAKE algorithm, which enabled use of a 2-fs integration time step.

TMA
A filament containing 52 peptides was constructed, and 5000 steps of ABNR minimization were computed. Before starting the simulation, one end was clamped in space by fixing the coordinates of the C carbons of four peptides at the end. Heating was performed in the same way as described above, followed by 400 ps of equilibration, during which temperature was rescaled in the 300 ± 10 K range. The production run was performed for 3 ns, and the coordinate trajectory was saved every 0.8 ps.

NMA
For a given filament structure (60 peptides in size), 200 steps of ABNR minimization were performed. The system was further relaxed by heating to 100 K in 40 ps and equilibrating for 100 ps, followed by full ABNR minimization. The diagonalization-in-mixed-basis, or DIMB (53) module, in CHARMM was applied to the minimized structure to calculate the normal modes. After 3000 iterations, the first 600 lowest-frequency modes were obtained.

SMD
A filament with a clamped end containing 44 peptides was constructed, heated, and equilibrated in the same way as in TMA. Three production runs, each lasting 4.4 ns, were performed during which a ramped external force was applied to the free end of the filament. In the first two cases, the forces were applied transverse to the filament axis, increasing from 0 pN to 150 pN, in 15 pN steps each lasting 400 ps. In the third simulation we applied a tensile force increasing from 0 pN to 300 pN with a step size of 30 pN.

Supramolecular packing geometry
We first constructed a series of ß-sheet tapes that are compatible with the filament dimensions obtained from AFM (48). Under AFM, the RAD16II filaments appear as straight, branched tapes 10 nm in width and 2 nm in height (Fig. 1 b). Considering the 3–5 nm radius of curvature of the AFM tip and the molecular dimensions from Fig. 1, we concluded that the filament is a ß-sheet bilayer (7), one molecule in width. Due to the alternating arrangement of hydrophobic and hydrophilic side chains, the ALA side chains were placed inside the bilayer (4).

Preliminary simulations suggested that the arrangement between peptides is antiparallel; the minimized energy of the parallel alignments was +15 kcal/mol per peptide higher than that of the ground state of the antiparallel ones—a significant difference, despite considering numerical accuracy. This is mainly due to favorable electrostatic interactions between charged side chains in antiparallel arrangements (48). We then considered all possible in-register antiparallel ß-sheet patterns. Due to the asymmetrical distribution of the backbone hydrogens and oxygens, there are two distinct antiparallel hydrogen-bonding patterns on each side of a peptide, which we call S1, S2, and S3, S4, respectively (Fig. 2). Alternating these patterns give antiparallel ß-sheets, where the sheet composed of Si (i = 1, 2) and Sj (j = 3, 4) is denoted as Sij, resulting in S13, S14, S23, and S24. For example, a ß-sheet composed of three peptides will be formed by adding a new peptide to the S1 pattern in Fig. 2 from below. In this case, the new peptide has the choice of making either S3 or S4 pattern with the dark peptide on the lower part of S1, resulting in S13 or S14 (48).

View larger version (20K): [in this window] [in a new window]   FIGURE 2  Four hydrogen-bonding patterns in a ß-sheet of two RAD16II. Hydrogen bonds are indicated by vertical dashes. Hydrophilic side chains are pointing out of the page. The distance between the peptides is 4.8 Å (see (48)).

View larger version (37K): [in this window] [in a new window]   FIGURE 3  Relative orientation between peptides in the ß-sheet bilayer. Each arrow represents a peptide, pointing from N- to C-terminus. Peptides in the upper layer are shaded.

View larger version (57K): [in this window] [in a new window]   FIGURE 4  (a) Ribbon description of the filament and its characteristic motion: (b) bend, (c) splay, (d) stretch, and (e) torsion. The coordinate origin is set at the center of the cross section at the left end of the filament, and oriented such that the z axis is along the filament axis (L), x/y axes are along the long/short (W/H) edges.

K

(1)

where and are the moments of inertia of the cross section with respect to the corresponding axis of deflection. Here the integrations are performed over the cross-sectional area of the filament, and x and y are distances measured along the short (long) axis passing through the center of the filament for splay (bend). For a rectangular cross section,

(2)

For an isotropic material, K_T and K are related to Y and through the relations

(3)

where A = WH is the cross-sectional area, Y/2(1 + ) is the shear modulus, and J is the polar moment of inertia (55). For a rectangular cross section,

(4)

with = W/2H. For torsion, the relevant moment of inertia of the cross section (see Eq. 6) is defined along the filament axis:

(5)

Validity of this linear elastic description will be examined a posteriori when we analyze the simulation results. In the sections below, we discuss specific simulations from which the above quantities can be determined.

Wave equations and dispersion relations
In TMA or NMA, the vibrational characteristics of the filament can be related to its mechanical stiffness. For small deformations, the displacement variables as functions of the axial coordinate z and time t for bend (u_B(z, t)), splay (u_S(z, t)), stretch (u_T(z, t)), and torsion (u(z, t)) satisfy the wave equations (54)

(6)

where _l (_v) is the mass per unit length (volume) of the filament. If there is a drag force caused by solvent medium, the left-hand side of Eq. 6 will involve first-order differentials (_tu and _t) (36). However, since our simulations were done either in zero viscosity continuum solvent (TMA) or in semi-vacuo (NMA), the wave description without dissipation is more applicable.

The general solution of Eq. 6 can be expressed as a linear combination of (hyperbolic) sinusoidal waves (56),

(7)

Dispersion relations between the wave number k and the angular frequency can be found by substituting Eq. 7 into Eq. 6:

(8)

Different boundary conditions for TMA and NMA determine specific linear combinations of the general solutions (Eq. 7). The wave number k and angular frequency are then limited only to discrete values, k_n and _n (n = 1, 2, ···), as shown below.

Vibration of a cantilevered filament
In TMA, the filament is thermally vibrating with one end clamped. The corresponding boundary conditions for bend or splay are (54): u_{B, S}(0) = u'_{B, S}(0) = 0 and u''_{B, S}(L) = u'_{B, S}(L) = 0, which give

(9)

Here the wave number (k_n) satisfies

(10)

The values of the first three modes are k_nL 1.8751 (n = 1), 4.6941 (n = 2), and 7.8548 (n = 3). In the absence of viscous drag, the amplitude a_{(B, S), n} of the n^th mode is, in principle, determined by the initial conformation of the filament. In simulations, however, the initial condition corresponds to thermalization of the filament, so that the filament only undergoes thermally induced motions.

For stretch and torsion, the boundary conditions are u_T,(0) = u'_T,(L) = 0, with the corresponding solution

(11)

In analyzing the simulation data, we traced the spatial coordinate and torsion angle of the free end (z = L) by the method explained in Measurement of Deformations. Time-domain Fourier transforms were performed using the FFTW package (57), which gave characteristic vibrational frequencies of individual modes. Combined with the wave numbers obtained from Eqs. 10 and 11, the values of all four constants of mechanical compliance (K_B, K_S, K_T, and K) can be determined via Eq. 8.

Normal modes of a freely vibrating filament
Near the ground state, the interaction potentials between atoms are assumed to be harmonic. Diagonalization of the corresponding harmonic interaction potential matrix (Hessian matrix) yields the normal modes of the system. Since the low-frequency modes represent the global motion of the system, we expect them to exhibit wavelike behaviors. Since it does not directly follow the motion of the filament as a function of time, NMA is a useful complement to TMA or SMD.

Unlike TMA or SMD, the filament is not clamped, so the corresponding boundary conditions are u''_B,S(0) = u'_B,S(0) = 0 and u''_B,S(L) = u'_B,S(L) = 0 (54), with the corresponding solution

(12)

Similar to TMA, the wave number (k_n) is given by the relation

(13)

with the values of the first three modes corresponding to k_nL 4.7300 (n = 1), 7.8532 (n = 2), and 10.9956 (n = 3). For stretch and torsion, the first spatial derivatives at both ends are zero (u'_T, (0) = u'_T, (L) = 0) so that

(14)

After calculating normal modes, we captured the maximally deformed filament conformations for individual modes, which occur at a quarter time-point of respective vibrational periods. Among these, we selected those corresponding to the four characteristic deformations. These allowed us to calculate the mechanical stiffness by use of Eqs. 8, 13, and 14.

Filament deformation by external force
Although TMA and NMA identify passive vibrations of the system, SMD more actively seeks its response by applying an external force. For a linear elastic rod, the displacement of the free end d and the applied force F satisfy the relations (54)

(15)

From the slopes of the force-displacement curves, the filament's compliance can be calculated. Also, linearity of the curve would be an a posteriori check for our initial assumption of the filament as a linear elastic material.

Measurement of deformations
Axial length and deformation
We identified the contour at discrete points u(z_i, t) (i = 1,2,3,···), by locally averaging the coordinates of C atoms on four successive peptides, two from each sheet. From this, the filament length L(t) is computed as

(16)

The average, L(t), was used to determine the wave number k_n in TMA and NMA. Fluctuation of L(t) was <0.2%, far smaller than other sources of errors, so we treated the filament length to be fixed, L = L(t), where ··· denotes time average.

Transverse and torsional deformations
In TMA and NMA, bend had the largest amplitude and could be measured simply by projecting u(z_i, t) onto e₂:

(17)

In SMD, the external force produced a much greater deflection than was observed in TMA or NMA, so we used the same approach as for splay: u_S(L, t) = u(L, t) · e₁. For TMA and NMA, splay and torsion had amplitudes <1/10 of those for bending. For a more accurate measurement, we introduce a local coordinate system {n_l(z)|l = 1, 2, 3} along the filament axis. We first set n₃(z_i) = (u(z_i+1) – u(z_i))/|u(z_i+1) – u(z_i)|, tangential to the local filament axis. We then chose n₁ to lie in the plane of the local cross section of the deformed filament, which fixes n₂ = n₃ x n₁, perpendicular to the ß-sheet. With these definitions, n_l(z_i+1) and n_l(z_i) are related by the linearized rotation matrix (56)

(18)

where _1,2,3 are functions of z_i. Among these, ₁ (₂) is related to local bend (splay), while ₃ represents local torsion. By solving Eq. 18, we get for splay, , and for torsion, . Typical values of _i are such that _i – sin(_i) < 0.02. Moreover, for both TMA and NMA, our analysis was based on vibrational frequency, not on amplitude, thus the error caused by the linearization in Eq. 18 is negligible.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX: BENDING STIFFNESS OF... ACKNOWLEDGEMENTS REFERENCES

 
Identification of the supramolecular structure
For each filament structure constructed via the procedure in Methods, the energy and ASA per peptide were measured (Fig. 5). In contrast to the previous case of another ß-sheet peptide filament (48), there was no distinct lowest-energy state. Although S1313 or S1324 were the lowest in energy, S2424 and S1423 had energy levels that fell within the range of numerical uncertainty. The ASA of these states were among the lowest as well, and all four maintain the initial straight configuration throughout MD, whereas less stable ones deformed (Fig. 6). Molecular width and height, measured as end-to-end distances, ranged between 4.80–6.23 nm (W), and 2.30–2.48 nm (H), respectively, comparable to those from the AFM image, W 6 nm, and H 2 nm (Fig. 1 b). As shown in the sections below, the elastic moduli of these filaments were also similar. Thus the exact register of hydrogen bonds does not seem to be important in determining single filament mechanics.

View larger version (23K): [in this window] [in a new window]   FIGURE 5  Relative stability of different filament structures. (a) Average conformational energy during MD, (b) minimized energy of the average configurations, and (c) solvent-accessible surface area. Values were measured on per peptide basis.

View larger version (25K): [in this window] [in a new window]   FIGURE 6  Comparison of the molecular configuration of different packing methods at the end of MD.

(19)

where m_i and r_i are the i^th atom's mass and distance from the filament axis. We averaged Eq. 19 over time in the case of TMA, and over vibrational period for NMA. Measured values of _vI were, for S1313, S1324, and S2424, 22.4–22.9 kDa·nm (TMA), 19.8–21.6 kDa·nm (NMA), and for S1423, 13.6 kDa·nm (TMA), 9.85 kDa·nm (NMA). Since _vI fluctuated <1%, this small error was ignored in the estimation of K.

TMA and NMA
Time-domain Fourier-transforms were used to analyze the deformations in TMA. Peaks in the frequency spectrum gave _n in Eqs. 9 and 11 (Fig. 7). Since the first mode (n = 1) exhibited the sharpest, most well-defined peaks, we used these for analysis (Fig. 7, insets). Even so, the peak width contributed far more to the uncertainty of the stiffness values than the variations in filament length or densities. The wave number k₁ was estimated from Eqs. 10 and 11.

View larger version (32K): [in this window] [in a new window]   FIGURE 7  Normalized frequency spectrum of the S1313 filament from TMA. (a) Bend, (b) splay, (c) stretch, and (d) torsion. Arrows indicate the first and the second modes. All other filaments show similar behaviors. (Insets) Closeup near the first peaks of the four filaments.

View larger version (14K): [in this window] [in a new window]   FIGURE 8  Frequencies of the 30 lowest normal modes. The first six modes were ignored, since they represent translation and rotation of the center of mass. The identified characteristic motions are marked by different symbols: bend (x), splay (), stretch (), and torsion ().

View larger version (28K): [in this window] [in a new window]   FIGURE 9  Vibrational conformations of the filament in NMA. (a) Bend, (b) splay, (c) stretch, and (d) torsion. The modes identified in Fig. 8 are sorted and plotted together. The mode number n (Eqs. 13 and 14) is in angular brackets. The axes are both normalized, with the x axis as the axial position and the y axis as the vibrational amplitude. Shaded lines denote the solutions to Eq. 12.

K

View larger version (21K): [in this window] [in a new window]   FIGURE 10  Stiffness of the filaments in TMA and NMA. (a) Bend, (b) splay, (c) stretch, and (d) torsion. The error bars were obtained from the width of the first peak in the frequency spectra (Fig. 7) in TMA, or by averaging different modes in NMA (Fig. 9). Since only one stretching mode was identified for each filament (NMA), there is no corresponding error bar in c.

View larger version (21K): [in this window] [in a new window]   FIGURE 11  SMD simulation of the S1313 filament (L = 95.1 Å; see Methods). (a) Bend, (b) splay, and (c) stretch. The error bar represents the root-mean-square fluctuation of the filament tip at each force level. Straight lines of a and b are linear fits, showing linearity of the response (Eq. 15). (Inset) SMD simulation of the AGA16 filament.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX: BENDING STIFFNESS OF... ACKNOWLEDGEMENTS REFERENCES

 
Supramolecular structure of the RAD16II filament
Our results indicate that RAD16II exhibits several similarly stable filament structures. Although this finding makes it difficult to predict with confidence the one most likely to occur, it also might help to explain the branches observed in AFM experiments (Fig. 1 b). As seen in Fig. 3, since S1313, S1324, and S2424 have a tilt angle different from that of S1423, a mismatch between these structures could lead to branching. To test this hypothesis, we tried different supramolecular packing geometries of another peptide with a similar sequence, RAD16I (RADARADARADARADA). Experimentally, RAD16I does not produce branches, but rather exhibits sharp bends or kinks (data not shown). Preliminary simulation showed that in the case of RAD16I, S1313 (–3106.78 ± 3.18 kcal/mol) and S1324 (–3106.01 ± 5.17 kcal/mol) are the dominant structures (the number in parentheses is the average conformational energy per peptide during MD, similar to those in Fig. 5 a). The structures with the next lowest energies are S2424 (–3096.74 ± 3.46 kcal/mol) and S1423 (–3098.56 ± 2.73 kcal/mol), clearly not as stable as the first two. Another peptide, KFE8 (FKFEFKFE), produces neither branches nor kinks, and exhibits only one dominantly stable structure in simulation (48). Branches or kinks could thus be generated by competing interactions between similarly stable structures (Fig. 12). Although it would be difficult to probe different packing patterns near branch points, such a possibility could be explored by considering relative stabilities between different patterns from which the frequency of branching may be predicted.

View larger version (75K): [in this window] [in a new window]   FIGURE 12  A possible branching geometry caused by mismatch between similarly stable filament patterns.

View larger version (26K): [in this window] [in a new window]   FIGURE 13  (a) Nonbonded interaction energy of charged side chains in RAD16I and RAD16II filaments. (b) Van der Waals energy of the nonpolar side chain (ALA) of RAD16II.

(20)

Using a method similar to that in TMA for calculating stiffness, variations in _B,S can be estimated to give the upper and lower bounds of . For NMA, all possible pairs between five bending and two splaying modes were averaged to determine the aspect ratio. The measured values of from TMA/NMA are 3.75 ± 0.40/4.47 ± 0.71 (S1313), 3.93 ± 0.62/4.45 ± 0.44 (S1324), 3.62 ± 0.84/4.48 ± 0.91 (S2424), and 3.00 ± 0.68/3.74 ± 0.42 (S1423), slightly larger than, but consistent with, the approximate ratio of 3 obtained from the AFM image (W 6 nm and H 2 nm). If we instead assume that the filament's Young's modulus in stretch direction is the same as that in either bend or splay, we get the unrealistic estimate of >> 3 or << 3. These possibilities were thus discarded.

To determine W and H individually, we need one more condition. From Eqs. 2 and 5 and using the relation _l/_v = WH, we obtain

(21)

Equations 19–21 then provided estimates for H/W from TMA/NMA (in units of nm): 1.59:5.94/1.25:5.58 (S1313), 1.51:5.95/1.25:5.57 (S1324), 1.64:5.64/1.29:5.76 (S2424), and 1.71:5.11/1.19:4.46 (S1423). Among these, S1423 was noticeably narrower than the others, although its height was similar, reflecting its slanted arrangement of peptides (Fig. 3). The overall averages excluding S1423 were H = 1.42 nm, and W = 5.74 nm. Although these are comparable to those from the AFM experiment, the thickness is somewhat narrow, considering the 1.3 nm height of a single peptide in a bilayer (Fig. 1). In Molecular Origin for the Continuum Elastic Behavior, we explain this based on the observation that the filament elasticity is mainly determined by the strong, short-ranged interaction between peptide backbones rather than the comparatively weak side-chain interactions.

Solvation effect and comparison between the three simulations
Values of K_B,S measured from TMA were consistently larger than those from NMA (Fig. 10, a and b). This can be partly due to the different solvent models used. Analytic continuum solvent (ACE2) at 300 K was used in TMA, while the distance-dependent dielectric constant model (RDIE) was used in NMA. Furthermore, since the NMA calculation was based on harmonic perturbation of the minimum energy conformation, no temperature dependence was incorporated except when measuring the amplitude of vibration after all the modes were found. As a control, TMA for S1313 was performed with RDIE instead of ACE2, which gave K_B = (0.79 ± 0.13) x 10^–26 Nm² and K_S = (1.36 ± 0.05) x 10^–25 Nm², closer to the NMA result. Since ACE2 can partly account for the hydration shell formed around the charged side chains, the filament is expected to have a larger resistance to bend or splay, which are accompanied by the change in hydrophilic area. On the other hand, there is no noticeable difference between K_T, in the two simulations (Fig. 10, c and d). Under stretch, the filament had far less axial deformation compared to other deformational modes. Torsion, by definition, does not change the contour length. Since the overall size of the hydrophilic faces does not change under torsion or stretch, the hydration effect of charged side chains is not as important as in bend or splay, making them less sensitive to the choice of solvent models. An explicit water simulation would capture the effects of fluid viscosity, which can be incorporated into the analysis by addition of damping terms in Eq. 6. However, as we have shown above, different solvent models affect the values of stiffness by up to a factor less than 2, thus our main results will not depend strongly on the choice of solvation models. As in many other polymer systems, the elasticity is mainly determined by the material properties of the filament itself, rather than the solvent. Unfortunately, explicit water simulation was not feasible for our system, which is composed of 1000 residues with a total simulation period longer than 10 ns.

In this implementation, although TMA and NMA use different solvation models, the same wave equation formalism is used for analysis. On the other hand, TMA and SMD share the same solvation models, while SMD was analyzed by using a simple beam equation description without multiple frequency modes. Each approach also has additional limitations. Data from TMA are the noisiest. As mentioned above, NMA is more reliable for higher modes, but it is based on harmonic perturbation of a minimized structure without explicit temperature dependence. Since SMD measures the averaged force-displacement relation, it suffers less from noise, but it requires the most extensive computing resource. Due to such shared features and distinct limitations, the three approaches are mutually complementary, together enabling a more reliable analysis compared to an approach based on a single type of simulation.

An alternative approach using the equipartition theorem
The wave description was the main formalism used here for analyzing filament motion. Alternatively, in the case of TMA, we can apply the equipartition theorem to analyze the motion of the free end (35,37,58). According to the equipartition theorem, the average energy of each vibrational mode in equilibrium is equal to k_bT/2, where k_b is Boltzmann's constant and T is temperature, thus

(22)

These give K_B = (1.86 ± 1.06) x 10^–26 Nm², K_S = (1.94 ± 0.75) x 10^–25 Nm², K_T = (1.06 ± 0.06) x 10^–7 N, and K = (0.99 ± 0.36) x 10^–26 Nm² (averaged over four filaments), consistent with those from other approaches except for K, which is approximately two-thirds of the value based on the wave equations. For the equipartition theorem to be valid, many cycles of vibrational motion must be monitored, whereas the wave description, in principle, needs only one cycle of vibrational motion. Thus our main analysis based on wave equations would be computationally more efficient.

Applicability to helical filaments
Many biofilaments, including amyloid fibrils and some self-assembled ß-sheet peptide filaments, exhibit helical geometry. Our present computational approach can be generalized to such cases. For a given helical filament, TMA or NMA analyses similar to those presented here can be used by subtracting the equilibrium curvature from the measured angles. For SMD, Kirchhoff's equation (59), which is a generalized description for an elastic rod of an arbitrary shape, can be numerically solved to calculate the force-displacement relation similar to Eq. 15.

Elastic properties of the filament
Linear elastic behavior
In SMD, computed force-displacement curves for bend and splay show linearity of the response within the range tested (Fig. 11, a and b). In NMA, we further used Eq. 8 to determine whether angular frequency () and wave number (k) satisfy the relation k² (bend and splay), and k (torsion). Plots of versus k in log-log scale gave k^1.67 for bend, k^1.74 for splay, and k^1.20 for torsion (Fig. 14). Considering the uncertainty in (Fig. 7), we find that the filaments behave approximately as a linear elastic material in bending, splaying, and torsional modes.

View larger version (12K): [in this window] [in a new window]   FIGURE 14  Dispersion relation (K) in NMA. (a) Bend, (b) splay, and (c) torsion. Here, (for bend and splay) and (for torsion) are dimensionless. The rescaling parameter 0 depends on the mode; 0 = 80.5 ns–1 (bend), 230 ns–1 (splay), and 97.7 ns–1 (torsion). All plots are in log-log scale.

View larger version (16K): [in this window] [in a new window]   FIGURE 15  Young's moduli measured respectively from bend and stretch. Error bars are based on the uncertainty in angular frequency (TMA) and in root-mean-square of different modes (NMA).

Dependence on system size
To test whether our results depend on the system size, we performed TMA for the S1313 filament composed of 40 peptides (20 peptides on each layer), which yielded elastic stiffness slightly smaller than those of the 52-peptide filament in our main analysis. The differences were 7% (bend or splay), 9% (stretch), and 18% (torsion). Smaller values of elastic stiffness possibly reflect the increased contribution from filament ends, which are likely to be more flexible. In fact, the cutoff length of nonbonded interactions in the simulation was 15 Å, the length of a ß-sheet approximately four peptides in size (each separated by 4.8 Å). The presence of the exposed end will thus affect at most up to four peptides from the end. Although such edge effects can cause errors in the estimated values of stiffness, the main conclusions of our work should not be strongly influenced by the system size.

Molecular origin for the continuum elastic behavior
Similar to surfactants, self-assembly of ß-sheet peptides is driven by its amphiphilic nature (60). Balance between attractive (hydrophobic) and repulsive (hydrophilic and steric) interactions determines the global aggregate morphology, while formation of backbone hydrogen bonds confers the paracrystalline order of the cross-ß structure. In Appendix, we theoretically investigated the possible contribution of amphiphilic interactions mediated by side chains to filament elasticity by ignoring the shorter-ranged backbone hydrogen bonds. Without hydrogen bonds, the peptide filament behaves like a surfactant bilayer except for its different molecular geometry. One can then use a simple representation of the amphiphilic interactions that was originally developed for surfactants or lipids (61). Let h be the backbone-to-backbone distance between the two ß-sheets in a filament, and be the surface tension of the hydrophobic side chains. In Appendix, we show

(23)

Most hydrocarbons have ranging between 20 and 50 mJ/m², with amphiphilic molecules having values in the lower end of this range (62). This approach has been successfully applied to predict the bending stiffness of a lipid bilayer (62). However, in our case, using W 6 nm and h 1 nm, K_B 10^–28 Nm², approximately two orders-of-magnitude smaller than the value obtained from simulations. Even if we use the upper bounds, h = 2 nm and = 50 mJ/m², K_B = 1.2 x 10^–27 Nm². Thus amphiphilic interactions alone cannot account for the observed bending stiffness of the filament. Although they drive the self-assembly, once the structure is formed, its elasticity is likely to be dominated by the backbone hydrogen bond network.

To further test this idea, we ran SMD on a model ß-sheet bilayer filament composed of the peptide AGA16 (AGAGAGAGAGAGAGAG). The AGA16 filament has the backbone and hydrophobic side chain configurations identical to those of the S1313 filament of RAD16II except that the charged side chains (ARG and ASP) are absent. Since it is only a model filament, we used RDIE to avoid possible instability caused by solvation. Its bending stiffness was measured to be K_B = (1.38 ± 0.03) x 10^–26 Nm², similar to that of the RAD16II filament (Fig. 11, inset). Together with the theoretical argument above, this result confirms that the backbone hydrogen-bonding network is mostly responsible for the filament's elasticity.

The above finding may account for the observed values of the thin cross-sectional geometry of the filament. Since the backbone interaction plays a dominant role in determining the elasticity of the filament, the cross-sectional geometry of the corresponding elastic ribbon will be determined mainly by the network of hydrogen bonds. Indeed, the backbone-to-backbone distance between the two sheets is h = 0.78 ± 0.05 nm (see H = 1.42 nm), and the average distance between the first H and the last O atoms on either side of the peptide backbone is 5.4 nm (see W 5.74 nm). Considering the presence of side chains and capping groups at the N- and C-termini, the height and the width of the corresponding ribbon representation are expected to be slightly larger than those defined by the backbone hydrogen-bond network.

Implications for larger scale behavior
Persistence length and buckling force
Persistence length l_p of a polymer chain is related to its flexural rigidity K_f (40): l_p = K_f/k_bT. In our case, since K_B is 1/10th of K_S, K_f K_B (0.5 – 2.0) x 10^–26 Nm², yielding l_p 1.2–4.8 µm. Also, the critical buckling force of a filament of length l with pivoted ends is given by F_c = ²K_f/l² (54). In a peptide hydrogel, l corresponds to the typical pore size, 100 nm, yielding F_c 10 pN. Since a typical cell adhesion site (like one integrin binding) can generate a force on the order of 1–10 pN, the RAD16II hydrogel will make a soft three-dimensional substrate; cells will be able to mechanically deform the network and navigate through the hydrogel without necessarily degrading it.

Relation to the macroscopic rheology
One of us (34) previously constructed a cubic strut model based on cellular solid theory (63), to interpret macroscopic rheology data. The estimated Young's modulus of a single strut assembled by an eight-residue peptide (EFK8) was 0.6–20 MPa. The analysis was based on the assumed strut thickness of 10–30 nm, considerably larger than the thickness of a single filament of EFK8, which is approximately one-half the size of RAD16II. However, we recently have found through electron microscopy that self-assembled peptide filaments including RAD16II can form bundles in high concentrations (unpublished). Therefore, the strut model may capture the behavior of the bundle rather than the individual filaments. The question is then whether the stiffness measured from rheology is dominated by fibers (either individual filaments or bundles), or by the entanglement effect. That the effective Young's modulus of a model filament calculated from the strut model is much smaller than that of the individual filaments measured here (Fig. 15) may suggest that the network elasticity is governed more by the entanglement effect. The strut model and the single filament model are complementary in the sense that the former probes the average effective contribution of the filaments or bundles in a network, while the latter directly deals with a single filament. It remains as a future work to develop a polymer dynamical model where the effective stiffness from the strut model can be calculated as a function of the peptide concentration and single filament properties such as persistence length, cross-linking, or bundling behavior.

Comparison with other biofilaments
Last,