In-Silico Investigations of Calcium Phosphate Mineralization in

Publication Date (Web): March 8, 2018 ... Our investigations provide insights into the role of EVs on calcium phosphate mineral nucleation and growth ...
2 downloads 3 Views 1MB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

B: Biophysical Chemistry and Biomolecules

In-Silico Investigations of Calcium Phosphate Mineralization in Extracellular Vesicles Rudramani Pokhrel, Bernard S Gerstman, Joshua D. Hutcheson, and Prem P. Chapagain J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b00169 • Publication Date (Web): 08 Mar 2018 Downloaded from http://pubs.acs.org on March 9, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

In-Silico Investigations of Calcium Phosphate Mineralization in Extracellular Vesicles Rudramani Pokhrel1, Bernard S. Gerstman1,2, Joshua D. Hutcheson2,3, Prem P. Chapagain1,2* 1

Department of Physics, 2Biomolecular Sciences Institute, and 3Department of Biomedical Engineering, Florida International University, Miami, FL 33199, United States

Abstract Calcification in bone, cartilage, and cardiovascular tissues involves the release of specialized extracellular vesicles (EVs) that promote mineral nucleation. The small size of the EVs, however, makes molecular level studies difficult, and consequently uncertainty exists on the role and function of these structures in directing mineralization. The lack of mechanistic understanding associated with the initiators of ectopic mineral deposition has severely hindered the development of potential therapeutic options. Here, we used multiscale molecular dynamics simulations to investigate the calcification within the EVs. Results show that Ca2+-HPO42- and phosphatidylserine (PS) complexes facilitate the early nucleation. Use of coarse-grained simulations allows investigations of Ca2+-PO43- nucleation and crystallization in the EVs. Systematic variation in the ion-to-water ratio shows that the crystallization and growth strongly depends on the enrichment of the ions and dehydration inside the EVs. Our investigations provide insights into the role of EVs on calcium phosphate mineral nucleation and growth in both physiological and pathological mineralization.

*Corresponding Author: Email: [email protected] Phone: +1 305 348 6266

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. INTRODUCTION Calcification is a bio-mineralization process that involves accumulation of calcium phosphate, hydroxyapatite and apatite crystals within extracellular vesicles (EVs)1. Physiological calcification initiates bone and tooth formation2, whereas pathological calcification can occur in soft tissues, such as in blood vessels and joint cartilage, leading to debilitating conditions including atherosclerosis, osteoarthritis, and nephrocalcinosis3-6. The mechanisms associated with mineral deposition and growth in bone tissues remain uncertain, and despite the fact that calcifying EVs (or matrix vesicles) were identified over 50 years ago7, questions remain on the role and function of these structures in the calcification process8. Similarly, the mechanisms of ectopic mineralization in cardiovascular tissues remain unclear9 even though cardiovascular calcification is the most significant predictor of cardiovascular morbidity10 and the leading cause of death in the world11. The lack of mechanistic understanding associated with the initiators of this pathologic mineral has severely hindered the development of potential therapeutic options. Therefore, insight into the mechanisms of calcium phosphate mineral nucleation and growth would fill knowledge gaps in both physiological and pathological mineralization12-13. It is known that the formation of apatite crystals in bone and teeth occurs through an amorphous calcium phosphate (ACP) precursor phase14-15, which is characterized by the formation of stable prenucleation clusters. The small size of these clusters presents challenges in understanding the structural features of biological calcium phosphate precursor aggregates. Recently, Habraken et al.16 elucidated structural details in various stages in the formation of ACP and found that calcium-mediated clustering of calcium triphosphate Ca2(HPO4)32- forms the basis of the amorphous phase, which eventually leads to octacalcium phosphate (OCP) upon continued calcium uptake. Using a model surface that mimics the biological environment for calcification, Dey et al.17 showed that surface-induced apatite crystals develop from the nucleation of ACP, formed by aggregation of prenucleation clusters. It has been shown both experimentally and by simulations that negatively charged groups in collagen fibrils can serve as the binding sites for Ca2+ and induce the oriented nucleation of apatite18-19. Calcium ion (Ca2+) influx is regulated by annexin ion channels and inorganic phosphate (PO43-) is released from the hydrolysis of PO43-containing substrates by the enzyme alkaline phosphatase (TNAP). A growing number of studies have shown that mineral deposition in pathological calcification follows similar pathways to that of physiological calcifications, and EVs play an 2 ACS Paragon Plus Environment

Page 2 of 20

Page 3 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

important role in each20-23. EVs serve as the primary sites for the first mineral deposits in a wide range of biological calcification processes, from intramembranous bone formation, developing dentin and deer antler to pathological calcifications in aortic valves and cancer

24

. Calcifying

EVs are 30-100nm lipid enveloped structures that form by direct budding at the plasma membrane or through intracellular trafficking mechanisms24-27 and can attach to extracellular matrix proteins such as collagen II, IV, and X28-29. Compared to their parent cells, calcifying EVs are rich in phosphatidylserine (PS) and sphingomyelin (SM) lipids and annexin proteins, all of which enhance the process of calcification22, 24. PS-containing membrane surface is considered to be a substrate for nucleation in EVs30-31 but exact mechanisms of PS-induced nucleation within EVs is still unclear. In this paper, we used molecular dynamics (MD) simulations at the atomistic as well as coarse-grained level to investigate calcification within the EVs. While all-atom simulations give important details of the formation of prenucleation complex and early aggregates, coarse-grained (CG) method allows simulations of much larger systems such as EVs for reasonable time scales. CG molecular dynamics (CGMD) has been successfully used to study various biological systems32-36 and for large systems, CGMD is very useful to investigate the dynamics for time scales in tens of microseconds. Morphological changes of vesicles such as fusion and fission have been successfully described using CGMD model35, 37. Similarly, CGMD model has also been used to study the osmotic imbalance in EVs38.

2. METHODS 2.1 All-atom Simulations. The force field parameters for the HPO42- and PO43- were obtained from CHARMM General Force-field (CGenFF)

39

. The Ca2+-HPO42- prenucleation

complexes were first obtained by simulating 10 Ca2+ and 9 HPO42- ions in solution. For all-atom simulations, a system with lipid bilayer and Ca2+- HPO42- ions, and another with lipid bilayer and Ca2+-PO43- ions were built with Charmm-GUI web-server40. The complex, asymmetric membrane contains various lipids (POPC: POPE: POPS: POPI: PSM: CHOL) in the following ratios: (39: 9: 5: 4: 23: 20) in the outer leaflet and (8: 36: 20: 10: 5: 21) in the inner leaflet, which are in accordance of the asymmetric complexity of the plasma membrane 41-42. We used a slightly enhanced ratio of POPS on the inner leaflet to mimic the rich POPS environment of EVs 24. 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The prenucleation complexes were placed below the lipid bilayer, and were arranged on six planes as shown in Fig. S1. With 36 prenucleation complexes in each plane, the system contained 216 prenucleation complexes (i.e. 216 Ca2+ and 648 HPO42- ions). Additional 491 Ca2+ ions were added to allow Ca2+--mediated association of prenucleation complexes, making a neutralized system with the final Ca2+/HPO42- ratio of ~1.1 in the solution. For the PO43- system, HPO42- was replaced with PO43- to allow possible formation of apatite crystals. Additional 815 free Ca2+ ions were added to neutralize the system as well as to increase the Ca2+/P ratio. The neutralized system contained 1031 Ca2+ ions and 648 PO43- ions, giving Ca2+/ PO43- ratio as ~1.6. The details of the lipid compositions as well as the ions used in all-atoms simulations of both HPO42- and PO43- systems are given in Table S1. The membrane and ion systems were solvated with TIP3 water model. Gromacs-5.1.1 (gpu-version)

43-44

with CHARMM36 force-field was

used to perform all-atom simulations. The solvated systems were first minimized for 10,000 steps and equilibrated using the six-step protocol as given by Charmm-GUI. All atom simulations were propagated with a time-step of 2 fs. A 12 Å cutoff was used for the Lennard-Jones and coulomb potential. Coulomb interactions were calculated using PME 45. The particles within specified distances were tracked using the verlet neighbor scheme with a straight cutoff. The LINCS algorithm was used to control the rings and stiff bonds. The Berendsen and the Parrinello–Rahman barostats were used to maintain the pressure at 1 bar during the equilibration and production runs respectively. A semiisotropic pressure coupling was used with a compressibility of 4.5× 10−4 bar−1. The temperature coupling was done using Nose-Hoover thermostat and the temperature was maintained at 303.15 K. The Visual Molecular Dynamics (VMD)

46

was used for trajectory visualization, analysis, and

rendering.

2.2. Coarse-grained Simulations. For enhanced sampling of the crystallization inside EVs, we set up Coarse-grained systems of EVs and Ca2+-PO43- ions. Coarse-grained simulations of EVs were performed with Martini force field in which approximately four heavy atoms are mapped onto a single interaction site (bead) 47-49. The coarse-grained model for the phosphate ion (PO43-) was generated using the type Qa bead with -3 charges in the Martini force-field. The EVs with Ca2+-PO43- ions were generated in a series of steps. First, an EV with a radius of 10-nm was built using Martini-maker plugin in the Charmm-Gui web server 50-52 with asymmetric lipid bilayer containing lipids POPC, POPE, DPSM, POPS, DPCE and CHOL in the ratio of 4 ACS Paragon Plus Environment

Page 4 of 20

Page 5 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(1:1:3:0:1:3) in the outer layer, (1:2:2:3:1:3) in the inner layer , and a 10Å radius of cylindrical water pore with restraints to stabilize the vesicle during equilibration. Since the lipid distribution in EV is slightly different than in the plasma membrane 24, we reduced POPC and POPE but enriched DPSM, CHOL and POPS. To encourage Ca2+-PO43--POPS complex, all POPS lipids were put in the inner layer. The detail of the lipid composition of the EV is given in Table S2. The system was solvated with non-polarizable water using Martini22 model. The system was neutralized using Na+ and Cl– ions at 0.15mM. The system was then minimized for 6,000 steps and equilibrated in six-steps for a total of 135 ns, followed by 20 ns of unrestrained production run. The pore restraint was decreased in each step of the equilibration, with the harmonic force constant of 0.5 kJ/mol/nm2 in the first equilibration step to 0 in the unrestrained production run. The resulting equilibrated EV system was stable, with the pore fully sealed. For creating the EV systems with the Ca2+-PO43- ions, only the EV structure (lipids bilayer) of the final frame was saved and Ca2+ and PO43- ions were then added inside the EV. The EV system with the smallest ion concentration contained 795 Ca2+ and 542 PO43- ions in the form of a lose cluster, giving Ca2+/PO43- ratio of ~1.47. Instead of just replacing the water beads by the ions, placing a lose cluster of Ca2+ and PO43- ions in the beginning minimizes the shrinkage in the interior volume due to the ion clustering and avoids possible vesicle collapse. The EV containing Ca2+-PO43- ions was neutralized with Na+ and Cl-, and solvated with the polarizable water model53-54 using the python script INSANE55, which is more accurate at representing water interactions in solutions. Additional advantage of using the polarizable water model is that it avoids unwanted water freezing. A 100-ns equilibration was performed, followed by a 1-µs of production run. Several EV systems with increasing ion concentrations were built from the final frame of the 1-µs simulation of the EV system described above by replacing the water beads with Ca2+ and PO43ions, keeping the Ca2+/ PO43- ratio close to 1.45. We define the ion concentration, i.e. the ion-towater ratio () within the EV, as follows:

=

              

(1)

With four water molecules in each water bead, the effective ion concentration is /4. The number of ions and water beads used in different EV systems are given in Table S3. Each EV system was energy minimized for 14,000 steps and equilibrated for 100 ns, followed by 1.2 µs of production run with Martini2.2. P force-field (for polarizable water model) using Gromacs-5.1.1(gpu-version)

43-44

.

Coarse-grained simulations were propagated with a time-step of 15 fs. The Lennard-Jones and coulomb potential were cutoff at 12 Å, and the coulomb interactions were calculated using PME45. 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Berendsen and the Parrinello–Rahman barostats were used to maintain the pressure at 1 bar during the equilibration and production runs respectively, with an isotropic pressure coupling with a compressibility of 4.5× 10−4 bar−1. The temperature coupling was done with a velocity rescale algorithm and the temperature was maintained at 303.15 K.

3. RESULTS AND DISCUSSION 3.1. Clustering of Ca2+-HPO42- and Ca2+-PO43- and the Role of Lipid Interactions. Calcium phosphate mineralization occurs in multiple stages16-17. First, prenucleation complexes loosely associate in prenucleation clusters, which can subsequently aggregate into ACP and associate with the membrane. Although ACP was generally considered to be composed of Ca9(PO4)6 56, recent studies have demonstrated the presence of HPO42- in the ACP57-58. Habraken et al.16 showed that the prenucleation clusters of calcium phosphate17 are composed of ionic association of prenucleation complexes [Ca(HPO4)3]4-. A prenucleation complex consists of a calcium ion bound with three hydrogen phosphate ions. We simulated Ca2+ and HPO42- ions in the solution, with a Ca2+/ HPO42- ratio of 1.1 and prenucleation complexes are found to form spontaneously in all-atom simulations. As shown in Fig. 1a, this prenucleation complex is consistent with the structure obtained from ab initio calculations16, 59. With enrichment of Ca2+ ions, the calcium triphosphate prenucleation complex forms the basis for the crystallization into octacalcium phosphate and apatite in intramembranous mineralization60. To investigate the membrane interactions of the prenucleation clusters, we set up a system with the plasma membrane and the prenucleation complexes as well as additional Ca2+ ions in solution, giving a Ca2+/HPO42- ratio of ~1.1. We performed a 100-ns all-atom simulation for this system as discussed in the Methods section. We found that these prenucleation complexes quickly aggregate to form several clusters. A typical Ca2+-HPO42- cluster is shown in Fig.1a. In each cluster, the ratio of Ca2+/HPO42- is ~1.0. The Ca2+/HPO42- clustering process from the all-atom simulation is shown in Movie S1 and summarized in Fig S1. To encourage OCP-like aggregates, we set up another system with plasma membrane and Ca2+ and PO43- in the solution with the Ca2+/ PO43- ratio of ~1.6. A 100-ns all-atom simulation for this system showed that Ca2+ and PO43- quickly form clusters with Ca2+/PO43- ratio close to 1.3 (Fig.1a). Both Ca2+HPO42- and Ca2+- PO43- clusters are found to aggregate on the membrane surface at the inner leaflet to form a densified cluster (Fig. 1b and Fig. S2). In both systems, a significant association 6 ACS Paragon Plus Environment

Page 6 of 20

Page 7 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

with the membrane surface is observed, with phosphate head-groups of PS binding with the Ca2+ ions to localize the cluster.

Figure 1. a) Complexes of Ca2+ with HPO42- or PO43- in solution: prenucleation complex (left), Ca2+-HPO42- cluster (middle), and Ca2+-PO43- cluster in the solution (right). In the prenucleation complex, one calcium ion (white sphere) is surrounded by three hydrogen phosphate ions (sticks). Similarly, in Ca2+- HPO42- and Ca2+-PO43- clusters, calcium atoms are represented as white spheres, phosphorous atoms as orange spheres, oxygen atoms as red sticks and the hydrogen atoms as white stick. b) Stabilization of the Ca2+-HPO42- clusters at the plasma membrane. Solvent is not shown for clarity. On the right, the POPS lipid interacting with a Ca2+HPO42- cluster is highlighted. To examine the role of each lipid type on the localization and stabilization of the cluster, we calculated the number of Ca2+ ions within 3.5 Å of lipids in the inner leaflet. Figure 2 shows that compared to other lipids, POPS clearly stands out as having significant affinity for binding Ca2+ in both systems. For the system with Ca2+- PO43-, each POPS lipid has contact with more 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

than one Ca2+ ion on average at the end of 100 ns. This is followed by POPE, which has contact with less than one Ca2+ ion on average. Despite the negatively charged nature of POPI, it is found to have the least affinity for Ca2+ binding. Due to fewer Ca2+ ions in the system for Ca2+HPO42-, the count per lipid is much smaller (Fig. 2a).

Figure 2. Average number of Ca2+ ions per lipids on the inner (lower) leaflet of the membrane bilayer with (a) HPO42- and (b) PO43- in the solution. 3.2. Coarse-grained Simulations of Ca2+-PO43- Crystallization in EVs. Since crystallization is a slow process and large EV systems are computationally demanding for allatom simulations, we used Coarse-grained Molecular Dynamics (CGMD) simulations with the Martini force field to investigate the crystallization in EVs. Coarse-grained (CG) simulations have been used to study various biological systems32-36 for microsecond timescales. CG simulations have provided important information on the dynamics and physical processes in EVs, including the morphological changes35,

37

and osmotic imbalance38. For our CG

simulations, we represented the Ca2+ ion and PO43- ion each by a single bead. The general Martini model of CG water is known to result in rapid, unwanted freezing of the CG water, with a freezing temperature close to 290 K48, 61. The presence of the lipid bilayer as well as the quick clustering of the ions in the solution had a compounding effect on water freezing, and the EV systems freeze rapidly. A way around the problematic freezing is to add antifreeze particles to the solvent, which are basically larger beads of CG water48. For our system, the problem persisted even with the addition of 10% anti-freeze particles in the solution. Therefore, we chose 8 ACS Paragon Plus Environment

Page 8 of 20

Page 9 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

a polarizable water model54 for the solvent. Unlike the general CG water model, the polarizable CG water bead contains orientational polarizability and allows a more realistic representation of the water interactions with charged particles54. For a slightly increased computational cost, the polarizable water model avoids the aphysical water freezing54 for T>285 K. With this model, we did not observe water freezing in any of our simulations. As seen in the general CG water model62, we note that the isotropic nature of interactions of single CG bead to represent the PO43ion can introduce inaccuracy in the nucleation/crystallization, as well as in the final crystal structure. Also, as in the general CG water, the nucleation and crystallization of Ca2+-PO43- is likely enhanced by the use of a single-bead representation of PO43- ion. We performed CG simulations for several EV systems with different concentrations of 2+

Ca

and PO43- ions. As shown in Fig. 3a, the CG simulation accurately reproduces the

prenucleation complexes observed in all-atom simulations (Fig. 1a). Most prenucleation complexes consist of three PO43- beads bound by one Ca2+ bead (Fig. 3a), consistent with the proposed model of the prenucleation complex16. The Ca2+-PO43- prenucleation complexes form and break regularly in a fluid-like fashion until a group of complexes start to cluster together. Prenucleation complexes can assemble together via free Ca2+ ions to form a globular nucleation cluster16. With three PO43- beads bound to one Ca2+ in the prenucleation complex, the Ca2+/ PO43ratio would be only 0.33. To encourage nucleation, as well as to bring the Ca2+/ PO43- ratio between 1.3-1.6, we included additional free Ca2+ ions in the solution. This allowed further clustering of Ca2+- PO43- that would eventually form a nucleation cluster. As discussed later, we found that the nucleation can only occur if the ion-to-water ratio is relatively high (ρ>0.25). Although the loosely clustered aggregates can start anywhere in the solution, their stabilization by POPS lipids seems to be important (Fig. 2a). POPS interaction with the cluster via Ca2+ ions provides immobilization of the cluster and enhances nucleation and propagation. Once the nucleation cluster is stabilized, propagation of crystallization is swift. This observation is consistent with experimental results63-64.

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 20

Figure 3. a) Prenucleation complexes in solution. One calcium ion is surrounded by three phosphate ions (left). Prenucleation complexes aggregate to form Ca2+- PO43- clusters in solution (middle), followed by the binding of the Ca2+- PO43- cluster with POPS lipids inside of the vesicles (right). The negatively charged head group of POPS binds with the calcium ions (white) to form the stable cluster, which gives the templets for crystallization. Following the nucleation (highlighted by a circle), crystallization can quickly propagate throughout the EV.

b) Color

coded mobility of the Ca2+ ions according to their rmsf to show nucleation (white circle on the left), propagation (middle) and full crystallization (right). Color codes: blue- high mobility, white- intermediate mobility and red- low mobility. The Phosphate ions are colored orange throughout. For clarity, water beads are not shown.

To locate the nucleation event and to understand the propagation of crystallization, we tracked the mobility of the Ca2+ ions by calculating their positional root-mean-squared fluctuations (rmsf). The ions are highly mobile in the solution, but their mobility is reduced significantly once they start to cluster or clump together. We calculated the rmsf for each of the 10 ACS Paragon Plus Environment

Page 11 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Ca2+ ions in the system for a sliding window of 6 ns. For visualization of the nucleation and propagation, we color-coded the Ca2+ ions, from blue to red, with blue representing the highest mobility and red representing the least. Figure 3b displays the formation of the nucleation cluster as indicated by the group of slow-moving Ca2+ (colored red), followed by propagation and crystallization. The process of crystallization is shown in Movie S2. Examination of the final crystal structure shows significant number of water beads interstitially trapped in the crystal, as shown in Fig. S3. The trapped water molecules may have implications in crystal growth65 but the fact that the CG water bead contains four water molecules and has similar size as that of the ions, the inclusion of water is likely exacerbated in the CG simulations.

3.3. Dependence of the Ion-to-water Ratio (ρ ρ) on Crystallization. The first EV-ion system we simulated consisted of Ca2+- PO43- ions and polarizable water within the EV in the ratio of ρ=0.2, which did not result in crystallization at 303 K for up to 1 µs of a CGMD simulation. To facilitate crystallization, we systematically increased the ion concentration (ρ) by replacing the water beads within the EVs with Ca2+ and PO43- ions while keeping the Ca2+/ PO43ratio the same. We found that for ρ>0.26, nucleation and rapid crystallization occurred at 303 K within 1 µs of CGMD simulations. Figure 4 shows the structures of the representative EV systems at the final frame of each simulation. The Ca2+ and PO43- ions are seen clearly ordered for ρ>0.26.

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Final structures of the EV systems with various ion concentrations. For clarity, water beads are not shown. To track the crystallization at different ion concentrations, we calculated the Lindemann Index (or Lindemann parameter)66, δ, which is defined as the relative root-mean-squared bondlength fluctuation and is given by

) 2 ! = ( #$# % 1' 1

,,-3,

〈+,-. 〉0 % 〈+,- 〉.0 〈+,- 〉0

where rij is the distance between the ith and jth ions and 4 50 is the time average.

12 ACS Paragon Plus Environment

Page 12 of 20

Page 13 of 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The Lindemann index is generally used to analyze the melting behavior of solid crystals. In a solid, the atoms have small vibrational amplitudes about their equilibration positions and therefore yield a small Lindemann index, whereas the enhanced mobility of atoms in a liquid gives a larger Lindemann’s index. The general criteria for a solid is δ 0.1, though it depends on the material67-68. For each EV system, we calculated the Lindmann Index for the Ca2+ ions for sliding windows of 6 ns for the MD trajectory. As shown in Fig. 5, δ is high and does not change with time during the entire simulation for systems with ρ0.26, δ shows a clear transition from an initial liquid phase to a solid-like phase as the crystallization propagates throughout the system. Interestingly, the nucleation event is found to occur in a seemingly random fashion with respect to time. Nucleation for ρ=0.27 occurred around 250 ns, whereas it took over 650 ns for nucleation to occur for ρ=0.29, with a similar timescale for ρ=0.266. Although no apparent temporal pattern for nucleation is observed in these trajectories, more detailed investigations of the first passage times obtained and averaged from multiple trajectories for each system may present trends. Fig. 5 shows that once nucleation starts, it propagates through the whole volume in less than 50 ns in our CGMD simulations.

Figure 5. Lindemann Index (δ) for different ion concentrations; δ >~ 0.15 for the liquid-like phase and δ < 0.02 for the crystalline structure. 13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Both the Lindemann index δ and the rmsf (Fig. S4) follow similar trends as they both measure the positional fluctuations (or the mobility) of the ions. Rmsf gives the average fluctuation of an ion from its own mean position for a given period of time, whereas the Lindemann index gives the fluctuation of ions relative to all other ions. Also, calculation of the energy shows that the interaction (electrostatic+VDW) between calcium and phosphate ions becomes much stronger with the increase in ion concentration. As displayed in Fig. S4, low ion concentrations (ρ