Large-scale Molecular Dynamics Simulations Research Articles

Simulated and experimental study of the chip deformation mechanisms of monocrystalline Cu

Monocrystalline Cu exhibits excellent electrical and signal-transmission properties due to its absence of grain boundaries, making it a critical material for the production of micro-machinery and micro-components; however, achieving ultrahigh precision and ultralow damage machining of functional devices using traditional techniques such as grinding and polishing is extremely challenging. Consequently, nanocutting has emerged as an efficient means to fabricate monocrystalline materials with complex surface characteristics and high surface integrity. Nevertheless, the macroscopic cutting theory of metal materials cannot be applied to nanocutting. Accordingly, in this paper, both simulations and experiments were conducted to examine the chip deformation mechanisms of monocrystalline Cu. First, large-scale molecular dynamics (MD) simulations were conducted to gain a comprehensive understanding of the deformation behavior during nanocutting. This included examining the influencing factors and the variation patterns of the chip deformation coefficient, cutting force, and minimum cutting thickness. Subsequently, nanocutting experiments were performed using a specially designed nanocutting platform with high-resolution online observation by scanning electron microscopy. The experimental results served to verify the accuracy and reliability of the MD modeling, as they exhibited excellent consistency with the simulated results. Although this work considered monocrystalline Cu, it is believed that the elucidated chip deformation mechanisms could also be applied to other face-centered-cubic metals. These results are of great value for advancing the understanding of the mechanisms of ultraprecision cutting.

Self-organizing maps of unbiased ligand-target binding pathways and kinetics.

The interpretation of ligand-target interactions at atomistic resolution is central to most efforts in computational drug discovery and optimization. However, the highly dynamic nature of protein targets, as well as possible induced fit effects, makes difficult to sample many interactions effectively with docking studies or even with large-scale molecular dynamics (MD) simulations. We propose a novel application of Self-Organizing Maps (SOMs) to address the sampling and dynamic mapping tasks, particularly in cases involving ligand flexibility and induced fit. The SOM approach offers a data-driven strategy to create a map of the interaction process and pathways based on unbiased MD. Furthermore, we show how the preliminary SOM mapping is complementary to kinetic analysis, with the employment of both network-based approaches and Markov state models. We demonstrate the method by comprehensively mapping a large dataset of 640 μs of unbiased trajectories sampling the recognition process between the phosphorylated YEEI peptide and its high-specificity target lck-SH2. The integration of SOM into unbiased simulation protocols significantly advances our understanding of the ligand binding mechanism. This approach serves as a potent tool for mapping intricate ligand-target interactions with unprecedented detail, thereby enhancing the characterization of kinetic properties crucial to drug design.

Stability and cavitation of nanobubble: Insights from large-scale atomistic molecular dynamics simulations.

We perform large-scale atomistic simulations of a system containing 12 × 106 atoms, comprising an oxygen gas-filled bubble immersed in water, to understand the stability and cavitation induced by ultrasound. First, we propose a method to construct a bubble/water system. For a given bubble radius, the pressure inside the bubble is estimated using the Young-Laplace equation. Then, this pressure is used as a reference for a constant temperature, constant pressure simulation of an oxygen system, enabling us to extract a sphere of oxygen gas and place it into a cavity within an equilibrated water box. This ensures that the Young-Laplace equation is satisfied and the bubble is stable in water. Second, this stable bubble is used for ultrasound-induced cavitation simulations. We demonstrate that under weak ultrasound excitation, the bubble undergoes stable cavitation, revealing various fluid velocity patterns, including the first-order velocity field and microstreaming. These fluid patterns emerge around the bubble on a nanometer scale within a few nanoseconds, a phenomenon challenging to observe experimentally. With stronger ultrasound intensities, the bubble expands significantly and then collapses violently. The gas core of the collapsed bubble, measuring 3-4nm, exhibits starfish shapes with temperatures around 1500K and pressures around 6000bar. The simulation results are compared with those from Rayleigh-Plesset equation modeling, showing good agreement. Our simulations provide insights into the stability and cavitation of nanosized bubbles.

Structure and transport properties of LiTFSI-based deep eutectic electrolytes from machine-learned interatomic potential simulations.

Deep eutectic solvents have recently gained significant attention as versatile and inexpensive materials with many desirable properties and a wide range of applications. In particular, their characteristics, similar to those of ionic liquids, make them a promising class of liquid electrolytes for electrochemical applications. In this study, we utilized a local equivariant neural network interatomic potential model to study a series of deep eutectic electrolytes based on lithium bis(trifluoromethanesulfonyl)imide (LiTFSI) using molecular dynamics (MD) simulations. The use of equivariant features combined with strict locality results in highly accurate, data-efficient, and scalable interatomic potentials, enabling large-scale MD simulations of these liquids with first-principles accuracy. Comparing the structure of the liquids to the reported results from classical force field (FF) simulations indicates that ion-ion interactions are not accurately characterized by FFs. Furthermore, close contacts between lithium ions, bridged by oxygen atoms of two amide molecules, are observed. The computed cationic transport numbers (t+) and the estimated ratios of Li+-amide lifetime (τLi-amide) to the amide's rotational relaxation time (τR), combined with the ionic conductivity trend, suggest a more structural Li+ transport mechanism in the LiTFSI:urea mixture through the exchange of amide molecules. However, a vehicular mechanism could have a larger contribution to Li+ ion transport in the LiTFSI:N-methylacetamide electrolyte. Moreover, comparable diffusivities of Li+ cation and TFSI- anion and a τLi-amide/τR close to unity indicate that vehicular and solvent-exchange mechanisms have rather equal contributions to Li+ ion transport in the LiTFSI:acetamide system.

Enhancing Photothermal Energy Transduction through Inter- and Intramolecular Interactions of Multiple Two-Photon Dyes Appended onto Calix[4]arene.

Organic dyes-based photothermal agents (OPTAs) have received increasing attention as alternative to inorganic materials due to their higher biocompatibility and extensive diversification. Maximizing nonradiative deexcitation channels is crucial to improve the photothermal conversion efficiency (PCE) of OPTAs. This is typically achieved through individual molecular design or collective enhancement using supramolecular strategies. Furthermore, photothermal therapy (PTT) generally relies on linear one-photon absorption of the light source by the OPTA, with less consideration given to nonlinear two-photon absorption (2PA) strategies, despite their potential benefits. Here, a synergistic strategy, which combines intramolecular and intermolecular quenching, is employed to maximize the photothermal efficiency of diphenylamino-substituted distyryl dicyanobenzene (DSB), an outstanding two-photon-absorbing chromophore. One to three DSB units have been introduced on the conic p-tert-butyl-calix[4]arene (CX), serving as a preorganizing platform to allow aggregate formation and promote intramolecular quenching within the multichromophoric systems. Importantly, the multichromophoric molecules had very high two-photon absorption capabilities with cross sections (δ2PA) reaching maximal values of 3290 GM at 810 nm. Experimental data accompanied by large-scale molecular dynamics simulations and time-dependent density functional theory calculations shed light onto the interaction mechanism in those multiple DSB-appended CX compounds to rationalize their optical properties. Then, the formulation with Pluronic F127 amphiphile yields water-dispersible nanoprecipitates (Nps), in which the PCE is further maximized and the photobleaching is reduced due to the combination of intra- and intermolecular quenching. The high two-photon absorption in the near-infrared (NIR) window associated with the high PCE of these nanosized OPTAs could serve as a basis to future in vivo 2P-PTT applications.

Intrinsic Ferroelectric Switching in Two-Dimensional α-In2Se3.

Two-dimensional (2D) ferroelectric semiconductors present opportunities for integrating ferroelectrics into high-density ultrathin nanoelectronics. Among the few synthesized 2D ferroelectrics, α-In2Se3, known for its electrically addressable vertical polarization, has attracted significant interest. However, the understanding of many fundamental characteristics of this material, such as the existence of spontaneous in-plane polarization and switching mechanisms, remains controversial, marked by conflicting experimental and theoretical results. Here, our combined experimental characterizations with piezoresponse force microscope and symmetry analysis conclusively dismiss previous claims of in-plane ferroelectricity in single-domain α-In2Se3. The processes of vertical polarization switching in monolayer α-In2Se3 are explored with deep-learning-assisted large-scale molecular dynamics simulations, revealing atomistic mechanisms fundamentally different from those of bulk ferroelectrics. Despite lacking in-plane effective polarization, 1D domain walls can be moved by both out-of-plane and in-plane fields, exhibiting avalanche dynamics characterized by abrupt, intermittent moving patterns. The propagating velocity at various temperatures, field orientations, and strengths can be statistically described with a universal creep equation, featuring a dynamical exponent of 2 that is distinct from all known values for elastic interfaces moving in disordered media. This work rectifies a long-held misunderstanding regarding the in-plane ferroelectricity of α-In2Se3, and the quantitative characterizations of domain wall velocity will hold broad implications for both the fundamental understanding and technological applications of 2D ferroelectrics.

Data-driven modeling of dislocation mobility from atomistics using physics-informed machine learning

Dislocation mobility, which dictates the response of dislocations to an applied stress, is a fundamental property of crystalline materials that governs the evolution of plastic deformation. Traditional approaches for deriving mobility laws rely on phenomenological models of the underlying physics, whose free parameters are in turn fitted to a small number of intuition-driven atomic scale simulations under varying conditions of temperature and stress. This tedious and time-consuming approach becomes particularly cumbersome for materials with complex dependencies on stress, temperature, and local environment, such as body-centered cubic crystals (BCC) metals and alloys. In this paper, we present a novel, uncertainty quantification-driven active learning paradigm for learning dislocation mobility laws from automated high-throughput large-scale molecular dynamics simulations, using Graph Neural Networks (GNN) with a physics-informed architecture. We demonstrate that this Physics-informed Graph Neural Network (PI-GNN) framework captures the underlying physics more accurately compared to existing phenomenological mobility laws in BCC metals.

Large-scale molecular dynamics simulation and aggregate behavior research on asphalt

In this comprehensive study, a molecular dynamics (MD) model consisting of approximately 150,000 atoms was developed within the Strategic Highway Research Program (SHRP) framework to investigate the behavior of AAA asphalt, focusing on asphaltene aggregate dynamics behavior over a simulation time span of 0.282 microseconds. The findings reveal that the model achieves a stable state in terms of both energy and microstructure without significant changes or the formation of densely clustered molecules, despite the extensive simulation period. Notably, the patterns of asphaltene aggregation were found to significantly influence the MD system energy, with variations ranging from loosely organized structures due to dispersed asphaltenes to more structured, energy-efficient clusters formed by encapsulated asphaltene pentamers. The size of asphaltene aggregates emerged as a pivotal factor, affecting their encapsulation and leading to distinct formation patterns. Furthermore, the research highlighted that increasing shear rates prompt the depolymerization of asphaltene aggregates, shifting from large-scale structures to smaller ones and adapting dynamically to the flow conditions by aligning parallel to the shear direction. This study underscores the critical impact of asphaltene aggregate behavior and shear rates on the structural integrity and distribution within the asphalt model, offering profound insights into the MD of asphalt materials.

Exploring the constitutive activation mechanism of the class A orphan GPR20.

GPR20, an orphan G protein-coupled receptor (GPCR), shows significant expression in intestinal tissue and represents a potential therapeutic target to treat gastrointestinal stromal tumors. GPR20 performs high constitutive activity when coupling with Gi. Despite the pharmacological importance of GPCR constitutive activation, determining the mechanism has long remained unclear. In this study, we explored the constitutive activation mechanism of GPR20 through large-scale unbiased molecular dynamics simulations. Our results unveil the allosteric nature of constitutively activated GPCR signal transduction involving extracellular and intracellular domains. Moreover, the constitutively active state of the GPR20 requires both the N-terminal cap and Gi protein. The N-terminal cap of GPR20 functions like an agonist and mediates long-range activated conformational shift. Together with the previous study, this study enhances our knowledge of the self-activation mechanism of the orphan receptor, facilitates the drug discovery efforts that target GPR20.

Molecular dynamics simulation for phase transition of CsPbI3 perovskite with the Buckingham potential.

The CsPbI3 perovskite is a promising candidate for photovoltaic applications, for which several critical phase transitions govern both its efficiency and stability. Large-scale molecular dynamics simulations are valuable in understanding the microscopic mechanisms of these transitions, in which the accuracy of the simulation heavily depends on the empirical potential. This study parameterizes two efficient and stable empirical potentials for the CsPbI3 perovskite. In these two empirical potentials, the short-ranged repulsive interaction is described by the Lennard-Jones model or the Buckingham model, while the long-ranged Coulomb interaction is summed by the damped shifted force method. Our molecular dynamics simulations show that these two empirical potentials accurately capture the γ ↔ β ↔ α and δ → α phase transitions for the CsPbI3 perovskite. Furthermore, they are up to two orders of magnitude more efficient than previous empirical models, owing to the high efficiency of the damped shifted force truncation treatment for the Coulomb interaction.

Molecular dynamics simulations of neutron induced collision cascades in Zr — Statistical modelling of irradiation damage and potential applications

Understanding irradiation damage involves a multi-scale and multi-physics approach, integrating data from experiments, simulations and phenomenological models. This paper focuses on its early stages, specifically neutron-induced collision cascades in zirconium, as nuclear-grade zirconium alloys are widely used in fuel assemblies. We have gathered and analysed a significant sample of results from high-fidelity, large-scale molecular dynamics (MD) simulations, employing existing interatomic potentials and the two-temperature model to account for electron–phonon coupling. Our data can be directly applied to higher-scale methods.Furthermore, we carried out a comprehensive statistical analysis of the features associated with the defect production, including the number of defects, their distribution and size of the affected area. As a result, we developed a generative model of collision cascades that is parametric, hierarchical and stochastic, i.e. it takes into account a statistical nature of the phenomenon, is interpretable and shareable. This model has been developed with three primary objectives: to provide a sufficient descriptor of a cascade, to interpolate data obtained from high-fidelity simulations, and to demonstrate that the statistical model can produce representative distributions of primary irradiation defects. The results can be used to generate synthetic inputs for models at longer length and time scales, to provide fast approximations that take into account the morphology of introduced defects, and in general to serve as a powerful analytical tool.

Accuracy, transferability, and computational efficiency of interatomic potentials for simulations of carbon under extreme conditions.

Large-scale atomistic molecular dynamics (MD) simulations provide an exceptional opportunity to advance the fundamental understanding of carbon under extreme conditions of high pressures and temperatures. However, the fidelity of these simulations depends heavily on the accuracy of classical interatomic potentials governing the dynamics of many-atom systems. This study critically assesses several popular empirical potentials for carbon, as well as machine learning interatomic potentials (MLIPs), in their ability to simulate a range of physical properties at high pressures and temperatures, including the diamond equation of state, its melting line, shock Hugoniot, uniaxial compressions, and the structure of liquid carbon. Empirical potentials fail to accurately predict the behavior of carbon under high pressure-temperature conditions. In contrast, MLIPs demonstrate quantum accuracy, with Spectral Neighbor Analysis Potential (SNAP) and atomic cluster expansion (ACE) being the most accurate in reproducing the density functional theory results. ACE displays remarkable transferability despite not being specifically trained for extreme conditions. Furthermore, ACE and SNAP exhibit superior computational performance on graphics processing unit-based systems in billion atom MD simulations, with SNAP emerging as the fastest. In addition to offering practical guidance in selecting an interatomic potential with a fine balance of accuracy, transferability, and computational efficiency, this work also highlights transformative opportunities for groundbreaking scientific discoveries facilitated by quantum-accurate MD simulations with MLIPs on emerging exascale supercomputers.

Fundamental influence of irreversible stress–strain properties in solids on the validity of the ramp loading method

The widely used quasi-isentropic ramp loading technique relies heavily on back-calculation methods that convert the measured free-surface velocity profiles to the stress–density states inside the compressed sample. Existing back-calculation methods are based on one-dimensional isentropic hydrodynamic equations, which assume a well-defined functional relationship P(ρ) between the longitudinal stress and density throughout the entire flow field. However, this kind of idealized stress–density relation does not hold in general, because of the complexities introduced by structural phase transitions and/or elastic–plastic response. How and to what extent these standard back-calculation methods may be affected by such inherent complexities is still an unsettled question. Here, we present a close examination using large-scale molecular dynamics (MD) simulations that include the detailed physics of the irreversibly compressed solid samples. We back-calculate the stress–density relation from the MD-simulated rear surface velocity profiles and compare it directly against the stress–density trajectories measured from the MD simulation itself. Deviations exist in the cases studied here, and these turn out to be related to the irreversibility between compression and release. Rarefaction and compression waves are observed to propagate with different sound velocities in some parts of the flow field, violating the basic assumption of isentropic hydrodynamic models and thus leading to systematic back-calculation errors. In particular, the step-like feature of the P(ρ) curve corresponding to phase transition may be completely missed owing to these errors. This kind of mismatch between inherent properties of matter and the basic assumptions of isentropic hydrodynamics has a fundamental influence on how the ramp loading method can be applied.

Learning dislocation dynamics mobility laws from large-scale MD simulations

By dispensing with all the atoms and only focusing on dislocation lines, the computational method of Discrete Dislocation Dynamics (DDD) gains greatly over Molecular Dynamics (MD) in simulation efficiency of metal plasticity. But whereas in MD dislocations follow natural dynamics of atomic motion, DDD must rely on a dislocation mobility function to prescribe how a dislocation line should respond to the driving force exerted on it. However, reflecting our still incomplete understanding of ways in which dislocations move, mobility functions presently employed in DDD simulations entail simplifications and approximations of limited or, worse still, unknown accuracy and applicability. Here we introduce a data-driven approach in which the dislocation mobility function is modeled as a graph neural network (GNN) trained on large-scale MD simulations of crystal plasticity. We apply our proposed approach to predicting plastic strength of body-centered-cubic (BCC) metal tungsten and show that, once implemented in a DDD model, our GNN dislocation mobility function accurately reproduces the challenging tension/compression asymmetry of plastic flow observed both in ground-truth MD simulations and in experiment. Furthermore, subsequently validated by MD simulations, the same function accurately predicts plastic response of tungsten under conditions not previously seen in training. By demonstrating its ability to learn relevant physics of dislocation motion, our DDD+ML approach opens a promising avenue to bringing fidelity of DDD models closer in line with direct MD simulations at a much reduced computational cost.

Atomic cluster expansion potential for large scale simulations of hydrocarbons under shock compression.

We present an Atomic Cluster Expansion (ACE) machine learned potential developed for high-fidelity atomistic simulations of hydrocarbons, targeting pressures and temperatures near and above supercritical fluid regimes for molecular fluids. A diverse set of stoichiometries were covered in training, including 1:0 (pure carbon), 1:4 (methane), and 1:1 (benzene), and rich bonding environments sampled at supercritical temperatures, hydrogen rich, reactive mixtures where metastable stoichiometries arise, including 1:2 (ethylene) and 1:3 (ethane). A high-fidelity training database was constructed by performing large-scale quantum molecular dynamic simulations [density functional theory (DFT) MD] of diamond, graphite, methane, and benzene. A novel approach to selecting structures from DFT MD is also presented, which allows for the rapid selection of unique DFT MD frames from complex trajectories. Comparisons to DFT and experimental data demonstrate that the presented ACE potential accurately reproduces isotherms, carbon melting curves, radial distribution functions, and shock Hugoniots for carbon and hydrocarbon systems for pressures up to 100GPa and temperatures up to 6000K for hydrocarbon systems and up to 9000K for pure carbon systems. This work delivers a potential that can be used for accurate, large-scale simulations of shocked hydrocarbons and demonstrates a methodology for fitting and validating machine learning interatomic potentials to complex molecular environments, which can be applied to energetic materials in future works.

Uncovering Nanoindention Behavior of Amorphous/Crystalline High-Entropy-Alloy Composites.

Amorphous/crystalline high-entropy-alloy (HEA) composites show great promise as structural materials due to their exceptional mechanical properties. However, there is still a lack of understanding of the dynamic nanoindentation response of HEA composites at the atomic scale. Here, the mechanical behavior of amorphous/crystalline HEA composites under nanoindentation is investigated through a large-scale molecular dynamics simulation and a dislocation-based strength model, in terms of the indentation force, microstructural evolution, stress distribution, shear strain distribution, and surface topography. The results show that the uneven distribution of elements within the crystal leads to a strong heterogeneity of the surface tension during elastic deformation. The severe mismatch of the amorphous/crystalline interface combined with the rapid accumulation of elastic deformation energy causes a significant number of dislocation-based plastic deformation behaviors. The presence of surrounding dislocations inhibits the free slip of dislocations below the indenter, while the amorphous layer prevents the movement or disappearance of dislocations towards the substrate. A thin amorphous layer leads to great indentation force, and causes inconsistent stacking and movement patterns of surface atoms, resulting in local bulges and depressions at the macroscopic level. The increasing thickness of the amorphous layer hinders the extension of shear bands towards the lower part of the substrate. These findings shed light on the mechanical properties of amorphous/crystalline HEA composites and offer insights for the design of high-performance materials.

Unified understanding of the impact of semiflexibility, concentration, and molecular weight on macromolecular-scale ring diffusion

Conformationally fluctuating, globally compact macromolecules such as polymeric rings, single-chain nanoparticles, microgels, and many-arm stars display complex dynamic behaviors due to their rich topological structure and intermolecular organization. Synthetic rings are hybrid objects with conformations that display both ideal random walk and compact globular features, which can serve as models of genomic DNA. To date, emphasis has been placed on the effect of ring molecular weight on their unusual behaviors. Here, we combine simulations and a microscopic force-level theory to build a unified understanding for how key aspects of ring dynamics depend on different tunable molecular properties including backbone rigidity, monomer concentration, degree of traditional entanglement, and molecular weight. Our large-scale molecular dynamics simulations of ring melts with very different backbone stiffnesses reveal unanticipated behaviors which agree well with our generalized theory. This includes a universal master curve for center-of-mass diffusion constants as a function of molecular weight scaled by a chemistry and thermodynamic state-dependent critical molecular weight that generalizes the concept of an entanglement cross-over for linear chains. The key physics is how backbone rigidity and monomer concentration induced changes of the entanglement length, interring packing, degree of interpenetration, and liquid compressibility slow down space-time dynamic-force correlations on macromolecular scales. A power law decay of the center-of-mass diffusion constant with inverse molecular weight squared is the first consequence, followed by an ultraslow activated hopping transport regime. Our results set the stage to address slow dynamics and kinetic arrest in different families of compact synthetic and biological polymeric systems.

Simultaneous improvement in the elastic and fracture properties of graphene-epoxy nanocomposites–a computational perspective

The incorporation of stiff nano-additives (such as graphene) into a relatively soft polymer material (such as epoxy) usually leads to an improvement in elastic properties (i.e., Young's modulus) at the expense of fracture properties (i.e., tensile strength and toughness). Despite over a decade of research in polymer nanocomposites, we still lack a clear understanding of their structure-property relationships, which limits us from enhancing both elastic and fracture properties concurrently. Here, we performed large-scale reactive molecular dynamics simulations to study the deformation and fracture of model graphene/epoxy systems under uniaxial tension. A computationally efficient reactive force field was developed for graphene-epoxy system, allowing covalent bond formation, and breaking, which is crucial to model cross-linking in model epoxy as well as fracture. It was found the mechanical properties of the nanocomposite are very sensitive to the strength of the graphene-epoxy interface. As expected, elastic modulus increases with the interfacial strength. However, there appears to be an optimal interfacial strength to enhance the tensile strength and toughness. This is due to stress concentration occurring near the graphene edges at high interfacial strength, which leads to premature fracture. We show that by appropriately selecting an intermediate interface strength, one can simultaneously improve the ultimate tensile strength, toughness and the Young's modulus of the nanocomposite epoxy at ultra-low (∼0.2 % by weight) loading fraction of graphene additives. Our findings highlight the critical importance of properly engineered additive-matrix interfacial strength to develop high-performing polymer nanocomposites that are concurrently stiff, strong as well as tough.

Role of dislocations on martensitic transformation temperatures and microstructure: A molecular dynamics study

Microstructure and defects strongly affect martensitic transformations in metallic alloys. Significant progress has been made in understanding the atomic-level processes that control the role of grain boundaries and precipitates in these solid-to-solid phase transformations. Yet, the role of dislocations and their structures on martensitic transformation temperature and the resulting microstructure remains unclear. Therefore, we used large-scale molecular dynamics simulations to study the forward and reverse transformation of a martensitic material modeled after Ni63Al37 under cyclic thermal loading. The simulations reveal that dislocations in the austenite phase act as one-dimensional seeds for the martensite phase, which is present at temperatures significantly above the martensite start value. We find a reduction in the dislocation density during cyclic thermal loading, which results in the increase in martensite and austenite transition temperatures, in agreement with experiments. Importantly, we extracted a critical martensitic nuclei size for developing stable domains and found that relatively low dislocation densities are needed to grow independent martensitic variants resulting in a multi-domain structure.

Heterogeneous Structure and Dynamics of Water in a Hydrated Collagen Microfibril.

Collagen type I is well-known for its outstanding mechanical properties which it inherits from its hierarchical structure. Collagen type I fibrils may be viewed as a heterogeneous material made of protein, macromolecules (such as glycosaminoglycans and proteoglycans) and water. Water content modulates the properties of these fibrils. Yet, the properties of water and the fine interactions of water with the protein constituent of these heterofibrils have only received limited attention. Here, we propose to model collagen type I fibrils as a hydrated structure made of tropocollagen molecules assembled in a microfibril crystal. We perform large-scale all-atom molecular dynamics simulations of the hydration of collagen fibrils beyond the onset of disassembly. We found that the structural and dynamic properties of water vary strongly with the level of hydration of the microfibril. More importantly, we found that the properties vary spatially within the 67 nm D-spacing periodic structure. Alteration of the structural and dynamical properties of the collagen microfibril occur first in the gap region. Overall, we identify that the change in the role of water molecules from glue to lubricant between tropocollagen molecules arises around 100% hydration while the microfibril begins to disassemble beyond 130% water content. Our findings are supported by a decrease in hydrogen bonding, recovery of bulk water properties and amorphization of the tropocollagen molecules packing. Our simulations reveal the structure and dynamics of hydrated collagen fibrils with unprecedented spatial resolution from physiological conditions to disassembly. Beyond the process of self-assembly and the emergence of mechanical properties of collagen type I fibrils, our results may also provide new insights into mineralization of collagen fibrils.