Molecular Dynamics Study of Collagen Fibrils: Relation between Mechanical Properties and Molecular Chirality

Collagen is a basic biopolymer usually found in animal bodies, but its mechanical property and behavior are not sufficiently understood so as to apply to effective regenerative medicine and so on. Since the collagen material is composed of many hierarchical structures from atomistic level to tissue or organ level, we need to well understand fundamental and atomistic mechanism of the collagen in mechanical response. First, we approach at exactly atomistic level by using all-atom modeling of tropocollagen (TC) molecule, which is a basic structural unit of the collagen. We perform molecular dynamics (MD) simulations concerning tensile loading of a single TC model. The main nature of elastic (often superelastic) behavior and the dependency on temperature and size are discussed. Then, to aim at coarse-graining of atomic configuration into some bundle structure of TC molecules (TC fibril), as a model of higher collagen structure, we construct a kind of mesoscopic model by adopting a simulation framework of beads-spring model which is ordina-rily used in polymer simulation. Tensile or compression simulation to the fibril model reveals that the dependency of yield or buckling limit on the number of TCs in the model. Also, we compare the models with various molecular orientations in winding process of initial spiral of TC. The results are ana-lyzed geometrically and it shows that characteristic orientational change of molecules increases or decreases depending on the direction and magnitude of longitudinal strain.


Introduction
So far, materials engineering and science have mainly focused on the development of strong and tough materials by using relatively hard and rigid substances, such as metal or ceramics [1] [2]. However, these substances are not always good to be integrated into human bodies. Therefore, exploration of soft and durable materials other than existing materials has been paid medical attentions and utility in surgery equipment has been required. Nowadays, in medical field, a strong and tough, but at the same time soft, material which will be fit to flexible human textures is being explored. Sometimes, it is considered better that these materials should disappear and be absorbed into original biological tissues during use. Therefore, it is nice thinking that the developing materials should be constructed by biologically the same constituent as human organs. From materials point of views, human body is composed of some materials categorized by two types of properties, that is, a hard or a soft tissue. In the category of the former, famous examples are bones and tendons (connecting a muscle to bone).
For example, tendons are very strong like bone, but at the same time they behave in flexible manner like muscle. We should focus on these types of human hard tissues, in order to research and to develop new types of artificial materials. It is well known that microscopically these hard tissues like bone are composed of small fibrous constituent, called "collagen" [3] [4]. It will be worthy to study the construction and mechanical properties collagen from very elemental scale since its characteristics surely come from the behavior in atomic scale [5].
Collagen is known as a fundamental substance in human and animal bodies.
In atomic scale, collagen is composed of narrow molecular chain called tropocollagen (TC), which is prospected as a novel and bio-originated element of materials. The structure of TC has been extensively studied theoretically and experimentally, but there has been little knowledge about mechanical behavior in connection with its relevant hierarchical structure. Therefore, in this study, by modeling structure of TC in atomic scale, we simulate numerically mechanical response of a single fiber of TC (TC fiber) and its bundles (TC fibrils). In particular, we will focus on mechanical flexibility, rigidity and shape-memory effect of TC fiber and fibrils. From experiment and theoretical approach, it is revealed that TC furnishes an unique triple helical structure (i.e. molecular chirality), where three long peptide chains (successively connected amino acid) are uniformly twisted each other and form fibrillar shape.
In this study, we will observe local and molecular twist angle of TCs when tensile or compressive deformation is applied to the structure and we will consider the origin of mechanical properties. So far, it has been a lot of researches concerning the deformation and strength of animal bones and other biological hard-tissue [6]- [11]. In those which have been researched, properties of the collagen have been characterized with the resolution at most μm size, by using finite element method [12] [13] and macroscopic particle method and one successfully predicts the mechanical response of such biological materials. Indeed, we can expect that those simulation results could be applied to a developing field of regenerative medicine and biomechanics (combination of biology and mechanics) application. However, essential property of collagen is surely based on much smaller size as of atomistic TC structure. From now on, we should link those macroscopic simulations to theoretical approach in atomic and molecular scale like our present study [14].
The point is that biological tissue or material is providing high performance such as ductility and high-strength using a multi-scale strategy. It is worth investigating to analyze the origin of stability and mechanical response of collagen from atomic scale to meso-or macroscopic scales by building up its hierarchical and multi-scale model, which is also relevant to most synthetic polymers. In considering the hierarchical structure of collagen fibrils [15], in this study, coarsegraining method is applied to the molecular dynamics (MD) simulation of a TC fiber. Aiming at multi-scale analysis, the bundle of TC fibers (Fibrillar TC: FTC) is also modeled, and then the size dependency of FTC for mechanical properties is investigated. In general, molecular orientation (MO) in constructing a single TC fiber greatly affects the function and strength of fibrillar materials made up of polymers [16]. So, in the present study, we will investigate how difference of MO angles of TC affects the behavior by comparing results obtained for different initial MO angle. It will reveal how mechanical response of TC fibrils depends on the triple helix of TC, that is, the unique chirality of its elements.
Before we construct coarse-graining MD model, we report on the research to full-atomic (all-atom) modeling of TC fiber and TC fibrils. It is because the reliability of multi-scale view certainly depends on the smallest scale reality, i.e. atomistic structure and behavior. The all-atom TC model is developed from open digital data of molecular structure, that is, protein data bank (PDB). We will conduct the simulation of all-atom model of single TC fiber by using versatile MD software (NAMD) [17] and discuss the response during tensile loading test in the depth of molecular scale, at first. This paper is organized as follows. First, the hierarchical structure of collagen from the stage of microscopic TC molecules is explained and the conceptual ways of numerical modeling are introduced along our study. Secondly, the way of all-atom modeling of TC molecule is shown. The atomistic modeling can express the chirality of peptide chains and their real structure. Then, to exemplify all-atom simulation behavior, tensile and compressive loading tests of the single TC are shown and the mechanical properties obtained there are discussed.
Then, as a main body to study collagen in this study, a single TC molecule and its bundle are modeled by using coarse-graining method. The structures of the multiple fibers are constructed and are analyzed. This leads to further discussion on the behavior of collagen tissues with many hierarchical stages in reality, by scrutinizing mechanical properties such as Young's modulus (tensile loading case) and buckling limit (compression case) depending on the number of components used those models. The interesting behavior associated with chirality change during deformation is clearly convinced. In particular, twisting and un-Journal of Biomaterials and Nanobiotechnology twisting behavior are found from the simulation and they are mathematically characterized by using general geometrical parameters such as torsion and curvature of the line segments in each peptide changes, by means of the Frenet-Selet's differential relations.

Hierarchical Structure of Collagen Fiber
It has been well known that collagen tissues which are found in animal bodies are one of fundamental biological and structural substances. Collagen constructs hard tissues of animal body such as bone, tendon and cornea. In viewing at microscopic region, it has a fibrillar or rod-shaped form [15]. Very different from the microstructure of other inorganic material such as metal, ceramics or even simple synthetic polymers, collagen is built up by many hierarchical structural stages, from a molecular scale to a biological organ scale. In each stage of size scales, fibrous elements of that size scale can be found.
It is known that the smallest structural unit of collagen has a molecular size, which is called tropocollagen (TC) [4]. The size of TC is regularly 1.5 nm in diameter and 300 nm in length. TC is a kind of polymer chain, i.e. peptide chain, which is made up of specific kinds of amino acid, e.g., glycine (Gly), Proline (Pro) and Hydroxyproline (Hyp). Basically, these amino acids are repeatedly polymerized many times into a single chain and finally they forms very long and massive chain molecule, which is often called polypeptide.
The interesting point of the structure of TC is that three polypeptide chains are smoothly aligned and holding each other in right-handed helical manner (i.e. triple helical structure is given). In experimental investigation using micrograph, it is also revealed that many bundles of TCs are aligned in their longitudinal direction having an offset of one fourth in the length of TC. It is also assumed that TCs are firmly bound together by cross-linked mechanism, where two neighbor TCs are attached each other near their end regions. Figure 1 shows the abstract of the hierarchical structure of TC, as explained here. Therefore, as crucial concept in multi-scale modeling of collagen, we should start from atomic scale, but the gap between each size scale is too large to build up straightforwardly. The larger scale directly relevant to organ size would be hard to model by the present atomistic modeling, just due to the limited degrees of freedom treatable in any cutting-edge computation equipment. So, it is also necessary to introduce some coarse-graining modeling which retains characteristics of atomistic-or peptide-level (molecular unit level) scale.

All-Atom Modeling of Tropocollagen (TC) and Exemplified Results
As shown above, collagen tissues are composed of TCs, which all have a fibrillar form. Each TC is constructed from atomic-scale peptide molecules. The purpose of this study is to understand the essential mechanism of mechanical response of TCs. Therefore, it is crucial to look into the structure of single TC from a microscopic scale to a mesoscopic scale which will include collective dynamics of TC fibrils. The best way to do this is that we build a physical model so as to observe the behavior of TC atom-by-atom, or possibly we go down to quantum mechanical level and conduct ab initio modeling and simulations. So far, atomistic scale simulation of TC molecules in the collagen tissues has been extensively conducted by many researchers. For example, Monti et al.  [22]. They also made mineralization model of collagen fibril [23] and made cross-linking models [24] to clarify strengthening or toughening mechanism of the collagen material, where they were using coarse-graining model to capture continuum-like mechanical properties. The cross-linking is one of crucial mechanism for the strength of collagen and the researchers including us are also focusing on it from microscopic point of view recently [25] [26] [27]. In fact, from the advent of MD study of the collagen then till now, during one or two decades, a lot of MD simulations have been conducted and continue applying in various conditions and adding modeling of relevant microstructural molecules [28]- [39]. Solely by atomistic modelling, it is hard to understand the collective nature of TC and collagen fibrils since the gigantic data produced by molecular dynamics or atomistic simulations are to be beyond our control. But, at the same time, such small-sized model should not be thrown away because it has accuracy on dynamical behavior of TC. In such context, anyway, we will try to model TC in all-atom method and to see how it works, first. Simulation of a single TC model is constructed based on all-atom modeling as explained as follows [40].
Molecular dynamics code we use here is NAMD (Nanoscale molecular dynamics) [17], in which many biological polymer or chemical substance can be simulated in straightforward way, once one decides what molecule will be dealt Journal of Biomaterials and Nanobiotechnology with. As an interatomic potential (in the field of bio-chemical molecular analysis, it is often called "force field"), we use CHARMM force-field [41] in which usual energetics of intra-molecular interactions such as for stretching, bending, torsion (including improper torsion) as well as inter-molecular interactions (van der Waarls and coulomb interactions). At this point, we should do confirm validity of the force-field, and we might rely on its long-standing improvement of data concerning structural stability and energetics of relevant molecules, which many researchers have been made so far. Of course, any force-field has not perfectly constructed until now, but we can expect that more sophisticated force-field will be available in the future. Since all-atom modeling now required is for confirmation of basic structure of TC molecules in relation to dynamics, the simulation configuration here is sufficiently satisfied for the purpose.
A molecular model of TC is shown in Figure 2. Length of real TC will be about 300 nm or such, but due to computational restriction for us, a short length model of TC, which is called "short-TC", is constructed at first. The atomic coordinate of "short-TC" model is strictly built from PDB (protein data bank [42]) structural data designated for TC, named "1CAG" [43]. The length of 1CAG structure is about 8 nm, in which 30-polymerized motifs of amino acid bases, Gly-Pro-Hyp, are contained. Usually the unit of TC should be Gly-X-Y [44], where X or Y can be an arbitrary amino acid but it is assumed that most X-Y will be "Pro-Hyp" [38]. But, since the version of CHARMM force-field we are using [45] does not include the function for Hyp base and, in fact, the sequence of "X-Y" is not always Pro and Hyp, we decide to apply "Pro-Pro" to "X-Y" and we consequently adopt a "Gly-Pro-Pro" motif as a repeated unit of three bases. This handling might alter the reaction of TCs, but existence of Pro and Hyp as a group of atoms differ not so much and we think qualitative discussion on dynamics will be possible.
The computational model of TC having actual length, where "Gly-Pro-Pro" units are repeating until they are polymerized 355 times, is called "Long-TC" model. The length of "Long-TC" model is almost 300 nm that is just the same as that of real collagen fibrils. We will compare results of the "Long-TC" model with those obtained in the "Short-TC" model. The diameters of those two TC models are equal and 1.5 nm. It means two models have different aspect ratio though they are both micro-sized fibers.
Tensile loading simulation is conducted by using steered (or interactive) molecular dynamics (SMD) method [46]. In the SMD method, we set some con-  Table 1.  the number of atomic interactions in "Long TC" are much larger than in "Short TC", so force per one interaction of the larger specimen come to small and it exhibits softer behavior. We can also find that the gradient of those curves depends on temperature. Figure 4 shows the temperature dependence of Young's moduli which are estimated by the gradient of stress-strain curves at the beginning stage (strain from 0 to 0.1). The reason for this temperature dependency will be explained that in higher temperature fluctuation of atomic motion is augmented, then interaction between atoms tend to be reduced in average, resulting in lower stiffness of the structures.  In concluding this section, we can confirm sufficient stability of single TC fibril and observe its sufficient resistance against the external tensile loading. Based on these microscopic and fundamental observations and understanding concerning the structure and the dynamics of TC molecule, we would like to move on to a larger-sized modeling, in which we can involve higher-level of structural hierarchy essential in collagen material.

The Coarse-Graining MD Model of Single Tropocollagen (TC) and Fibril of TCs
In the context of multi-scale simulation of materials, some coarse-graining methods of atomic simulations have been discussed extensively so far. For example, coarse graining molecular dynamics (CGMD) [47], dissipated particle dynamics (DPD) [48], Particle modeling [49] have been already invented. As for collagen also, there have been some attempts of coarse-graining MD modeling [19]. It is indispensable that any coarse-graining method includes phenomenological view point. As the first step of modeling, focusing on reproduction of the characteristics of the material is effective and crucial. We use a simple method of grouping atomic or molecular masses, retaining their microscopic nature, i.e., outline conformation of molecular chains as a whole, inter-and intra-molecular interaction such as stretching, bending and so on. The framework of the present coarsegraining utilized to TC molecules here is often called "beads and spring" model [50]. It has been utilized for various polymer chains and we find it has been sufficiently successful in simulating polymer materials. In this method, successive six units of peptide chain, each composed of three bases, Gly, Pro and Hyp in order (in all-atom model, Hyp was technically replaced by Pro, though), are grouped into one large massive bead (all relevant masses are concentrated onto the center of one particle). They are connected each other by inter-particle potential functions (like springs between masses). The potential functions are of intra-molecular and inter-molecular parts similar to all-atom modeling. The intra-molecular potential function represents the interaction inside a molecule and comprises two-body (corresponding to stretch and contraction) and three-body (corresponding to bending) functions. The inter-molecular potential function represents cross-linking (CL) between ends of collagen fibers, as well as ordinary inter-molecular interaction (those by Van der Waals or electrostatic bonding). For inter-molecular interactions, conventional Lennard-Jones (LJ) type function form is adopted. These potential functions are shown in Equations (1)-(4). The potential parameters of TC molecule we use here were suggested in the literature [22].
K In the formulation of LJ-type potential, as usual, ε indicates a binding energy and σ corresponds to inter-particular equilibrium distance when interaction force changes from the repulsive to the attractive. In the intra-molecular functions, str k and bend k work as spring constants for the change of elongation ( ) 0 r r − and three-body angle ijk θ , respectively, where 0 r and 0 θ are equilibrium two-body distance and three-body angle, respectively. In the CL potential, β is an arbitrary constant for cross-linking (here we use β = 12.5 referring former studies [51]. In these coarse-graining modeling, we should exclude torsion potential. In molecular scale, torsion potential of polymer might be considered among four successive atoms, but, in coarsened scale, the effect of molecular torsion on larger chains will be not crucial. Besides, we suppose that torsional rigidity when it is twisted could be compensated for by geometrically nonlinear nature of stretching, bending and inter-particle potential we are including in this model. The actual method to build up the structure of TC and its fibril is shown in   In order to discuss the effect of molecular orientation (MO) on properties, we modify the basic model by changing the pitch length in making helix of each TCs. All those models have different angles of chains with regard to longitudinal (loading) direction of the fiber. Those angles can be called orientation of molecules as for fibril direction, so we name those models orientation-angled fiber of tropocollagen (OFTC) models. Those models are specified by the single parameter, orientation angle φ, as shown in Figure 6. Figure 6 also shows how to make three-dimensional triple helical structure of TC from the planar development picture of inclined alignment of three polypeptide chains. The geometric parameters of the computation model and calculations conditions are summarized in Table 2.   clearly exhibit different behavior. The specimen seem to cause totally twist of TC structure. Figure 8 shows Young's moduli which are estimated from stress-strain relation up to nominal strain of 10%. As the number of TCs increases, Young's modulus increases because the interaction point used by cross-linking particles will increase. With regard to temperature effect, Young's modulus becomes lower in higher temperature since for those nano-sized structures thermal fluctuation due to large kinetic energy directly enhance average distance between particles. As a result, tensile rigidity as the whole structure becomes lower in high temperature. Note that, in the present modeling, the temperature is for every beads calculated and is not equivalent to that for all-atomic degree of freedoms.

Mechanical Response in Tensile and Compressive Loading
But, since it corresponds to kinetic energy of total system, it can be used as an extent of thermal energy of TC fibrils or collagen material has.
In compressive loading, stress-strain relation shows a sudden drop. It is assumed that buckling event occurs at this point. Figure 9 shows the comparison between the number of TCs and the stress value in longitudinal direction when buckling occurs. The larger the size of TCs bundle is the higher buckling limit it has. Mechanically speaking, the larger size in diameter means increase of the second moment of the area against bending deformation it causes theoretically large buckling load (i.e., Euler's theory of elastic buckling). Similar to tensile loading, in higher temperature, buckling stresses become lower due to the decrease of intermolecular binding energy in average.
Let us consider the effect of the orientation of molecules inside TC fibril on loading behavior, by using OFTC models which are built with different molecular orientation (MO) angles inside TC fibril. Figure 10 and Figure 11 show stress-strain curves obtained by tensile or compressive tests, respectively, in temperature of 100 K for OFTC models. In tensile loading, a certain drop of stress occurs at around strain range from 0.1 to 0.4 (marked by circles in Figure 10).
The state when the stress shows such drop is recognized as yielding stress for the total specimen. The smaller MO angle φ is, the larger the yielding stress becomes. It means that, when three peptide chains in each TC molecule are wound loosely, the resistance against stretching becomes higher.    Likewise tensile case, the smaller the MO angle φ is, the larger the buckling stress become. In TC fibers, when the MO angle is small, all the molecules tend to align to longitudinal direction. In this situation, the directions of covalent bonding which will produce the largest bonding energy among interactions prone to coincide with tensile or compressive axis. Consequently, for smaller φ, the resistance against axial loading becomes larger.

Geometric Analysis of the Conformation of TC Molecules (as for Change of Molecular Orientation Angle)
Chirality can be defined as twist angle along a TC molecule. In OFTC models, each of line segment connecting neighbor beads is continuously identified along respective peptide chains. When three peptide chains around the outer surface of that cylinder is recognized, the angle with respect to the center line of the cylinder is regarded as the orientation angle. By calculating these angles in three dimension space, we will be able to evaluate the degree of local twist of molecular chain inside TC structure. We try to use a mathematical handling by differential geometry as follows. A set of equations shown in Equation (5) (5) where the vector T is tangential to the connecting vector between neighboring particles (unit tangent vector). The vector N is always perpendicular to T (unit principal normal vector), and the vector B (unit binormal vector) is perpendicular to both T and N. And scalars κ and τ are called curvature and torsion, respectively. These scalar factors, κ and τ, are calculated from positions of particles inside chiral chains in our MD simulations. We can propose that they quantitatively estimate the degree of chirality (twisting) and its change for respective chains.
In particular, value of torsion τ indicates a degree of twist of the three-dimensional curve. We will monitor the change of τ, or alternatively (1/τ) which means torsion radius, during deformation. Figure 12(a) and Figure 12 τ are not equivalent because molecular orientation (MO) angles in making each models are different. Moreover, these radii also alter during deformation. Overall trend is that torsion radius increases for tensile loading and decreases for compressive loading. It is understood that tensile deformation of the specimen helps all molecules to align in parallel direction Journal of Biomaterials and Nanobiotechnology to the tensile axis and the torsion value τ is prone to decrease during that aligning. On the other hand, the torsion value τ increases during compressive deformation. That is because positions of beads tend to move in circumferential direction rather than in longitudinal direction during compression. This mechanism is also valid for buckling in compression though it occurs suddenly. As shown in strain range 0~0.05 ε = ± in the graphs of Figure 12, the magnitude of torsion increment 1/τ is generally larger for compression rather than for tensile case. In tensile case, constraint from neighbor particles limit their move, but in compression particles are relatively allowed local moves and chains deforms up to relatively large amount. It is guessed that such local move noted here is, for example, slip motion between TC molecules.

Conclusion
In this study, mechanical properties of collagen are studied from the simulation of a microscopic and biological polymer material called tropocollagen (TC) by using atomic simulations constructed with multi-scale point of view. The simulations are two-fold: the first one is constructed by using all-atom modeling and molecular dynamics (MD) method, where we observed a precise behavior and mechanical properties of a single TC molecule having a unique triple-helical structure made up of polypeptide chains. Secondly, as main simulation method, a TC molecule and its fibril are modeled by adopting coarse-grained framework to MD simulations, retaining chain and cross-linking structures by using beadsspring model and their force field. Elastic moduli and buckling limits are compared as for the size of the fibril and for temperature conditions. Besides, it is recognized that molecular chirality (geometrically evaluated as torsion around the axis of tensile or compressive loading) of each peptide chains related closely to mechanical state of total TC fibril. If a coarse-graining method can be directly built upon the results obtained all-atom modeling and can be connected larger continuum-level analysis like finite element analysis (FEA), such simulation will clarify the total mechanics of collagen material in all sizes. We could make the basis for that.