Machine-Learning Enabled New Insights into the Coil-to-Globule

Subscriber access provided by University of Sunderland

Chemical and Dynamical Processes in Solution; Polymers, Glasses, and Soft Matter

Machine-Learning Enabled New Insights into the Coil-to-Globule Transition of Thermosensitive Polymers Using a Coarse-Grained Model Karteek K. Bejagam, Yaxin An, Samrendra Kumar Singh, and Sanket A Deshmukh J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.8b02956 • Publication Date (Web): 29 Oct 2018 Downloaded from http://pubs.acs.org on October 31, 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 Letters

Machine-Learning Enabled New Insights into the Coil-to-Globule Transition of Thermosensitive Polymers using a Coarse-Grained Model Karteek K. Bejagama, Yaxin Ana, Samrendra Singhb, Sanket A. Deshmukha* a

Department of Chemical Engineering, Virginia Tech, Blacksburg, Virginia 24061, United States b

CNH Industrial, Burr Ridge, Illinois 60527, United States

AUTHOR INFORMATION Corresponding Author * [email protected], Phone: +1 540-231-8785 Abstract We present a computational framework that integrates coarse-grained (CG) molecular dynamics (MD) simulations and a data-driven machine-learning (ML) method to gain insights on the conformations of polymers in solutions. We employ this framework to study conformational transition of a model thermosensitive polymer, poly(N-isopropylacrylamide) (PNIPAM). Here, we have developed first of its kind, a temperature-independent CG model of PNIPAM that can accurately predict its experimental lower critical solution temperature (LCST) while retaining the tacticity in the presence of explicit water model. The CG model was extensively validated by

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry Letters 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 2 of 20

performing CG MD simulations with different initial conformations, varying the radius of gyration of chain, the chain length, and the angle between the adjacent monomers of the initial configuration of PNIPAM (total simulation time = 90 µs). Moreover, for the first time, we utilize the non-metric multidimensional scaling (NMDS) method, a data-driven ML approach, to gain further insights into the mechanisms and pathways of this coil-to-globule transition by analyzing CG MD simulation trajectories. NMDS analysis provides entirely new insights and shows multiple metastable states of PNIPAM during its coil-to-globule transition above the LCST.

TOC GRAPHICS

The coarse-grained (CG) molecular dynamics (MD) simulations have emerged as powerful technique to model polymers, biomolecules, and various architectures of these materials.1–8 CG models that represent a group of atoms as beads allow the use of higher timestep (5 to 50 fs) as compared to the all-atom (AA) MD simulations (1 to 2 fs). This enables us to reach the time- and length-scales of several 10s of microseconds (µs) and 10s of nanometer (nm), respectively.9–12 Various methods like force-matching, iterative Boltzmann approach etc. have been utilized to develop CG models of these macromolecules and solvents.11,13–17 Recently, we have developed a novel approach that integrates artificial neural network (ANN) based machine-learning model and particle swarm optimization (PSO) algorithm with MD simulations to accelerate the development of CG models (ANN-assisted PSO method).18 Despite of all these advancements, developing CG models of any stimuli-sensitive polymers that can accurately capture their conformational transitions has been a grand challenge. This is mainly because capturing the delicate changes in the inter- and intra-molecular interactions between polymer and solvent, in response to change in the surrounding environment, is very difficult.19 A CG model that can capture these interactions,


2

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


if developed, can be used to study the architectures and self-assembly processes in stimulisensitive polymers. Indeed, there have been a few recent efforts made to develop CG models of thermosensitive polymers such as Poly(N-isopropylacrylamide), poly(ethylene glycol) etc.20–22 In recent years, a range of different machine-learning (ML) methods such as deep learning, manifold learning etc. have been utilized to gain insights on crystal structures, folding pathways for biomolecules, and aggregation of small organic molecules.23–26 A non-linear manifold learning is related to the algorithmic techniques of dimensionality reduction and it consists of different methods including diffusion maps, non-metric multidimensional scaling (NMDS) method, isomaps, Hessian Eigen maps etc.27 The NMDS is a data-driven ML scheme and it has been employed to reduce the high-dimensions in the MD trajectories of protein folding while preserving the crucial features.25 Use of such manifold learning method to analyze the CG MD simulation trajectories of stimuli-sensitive polymers may allow one to capture the rare events, if any, occurring during their conformational transitions. Thus, it can provide insights on the mechanism and pathways of conformational transition. Poly(N-isopropylacrylamide) (PNIPAM), a model thermosensitive polymer, exhibits a coil-to-globule conformational transition above its lower critical solution temperature (LCST) of 305 K.28 As this coil-to-globule transition alters its physicochemical properties, it finds applications in the field of biomedicine, energy, coating etc.29–31 Developing a CG model of PNIPAM that accurately predicts its LCST is very challenging, as the model should precisely capture both the polymer-polymer and polymer-water interactions. To the best of our knowledge, only three CG models of PNIPAM have been reported in literature.32–34 However, they have one or more of the following drawbacks: (i) Although they show a coil-like and a globule-like state below and above the LCST of PNIPAM, respectively. The exact LCST of these models is not reported.32 (ii) They are temperature-dependent models i.e. a set of FF parameters used to define the interactions between polymer and water, below and above the LCST are different.33 (iii) These models could not retain the chain tacticity. Note, the tacticity of PNIPAM is known to play an important role in determining the LCST of PNIPAM.34 and (iv) The CG MD simulations utilized the implicit water model.34 Hence, none of these models can capture the polymer-water and polymer-polymer interactions accurately.19


3


Page 4 of 20

Here, we present a novel computational framework that integrates CG MD simulations and data-driven machine-learning methods to gain insights into conformations of polymers in solution. We employ this framework to study the coil-to-globule transition in PNIPAM. We began by developing first of its kind temperature-independent, chemically and structurally accurate, transferable CG model of PNIPAM by coupling MD simulations with the PSO algorithm.18,35–37 We show that the PNIPAM CG model retains its tacticity and is able to reproduce its experimentally reported LCST behavior. Then, we employ NMDS method to analyze the CG MD simulation trajectories and to explore the mechanism and pathways of coil-to-globule transition in PNIPAM.25 The NMDS analysis shows multiple metastable states of PNIPAM, which has not been reported for any thermosensitive polymer. The first step of CG model development of PNIPAM involved the force-field (FF) optimization of its monomer, N-isoprolylacrylamide (NIPAM). The FF parameters were optimized to reproduce the properties obtained from the AA MD simulations. The development of NIPAM CG model was based on the features and properties obtained from the AA NIPAM simulations that are modelled using CHARMM FF parameters.38 The AA MD simulations of NIPAM molecules was carried out at 300 K (See Supplementary Information Section 1.1). The CG NIPAM model was mapped on its functional groups, which resulted in a model with three beads, namely C2M, AMP, and ISP that encompass the CH2-CH2, H-N-C=O, and isopropyl groups of a NIPAM monomer (Figure 1a). Beads connected with bonds interact with harmonic bond potentials (See Supplementary Information Section 1.2). Beads separated by two bonds interact via harmonic angle potential as well as through nonbonded interactions. All three CG beads in NIPAM are chargeless and interact with each other via 12-6 Lennard-Jones (LJ) nonbonded potential. The AA simulation trajectory was converted into CG trajectory based on the mapping scheme illustrated in Figure 1a. The equilibrium bonds and angle between different CG beads were determined based on the mean value of the distributions obtained from the mapped AA trajectories. The values of σ for AMP and ISP beads were chosen based on the radial distribution functions (RDF) obtained from mapped AA trajectory. Nonbonded interaction parameters of C2M were adopted from our recently developed CG hydrocarbon model.36 Thus, to develop a CG model of NIPAM, we had to optimize the following five FF parameters: kb (for C2M-AMP; AMP-ISP), kθ (C2M-AMP-ISP), εAMP, and εISP. These FF parameters were optimized to reproduce the bonded distributions, RDFs obtained from the mapped trajectory and its properties such as density, surface tension, and heat


4

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


of vaporization obtained from NIPAM AA MD simulations. The optimization of FF parameters was performed by coupling PSO method with MD simulations. The timestep of 15 fs was used to perform the CG MD simulations at 300 K by using NAMD simulation package.39 The optimized FF parameters for NIPAM CG model are listed in Table S1 of the Supporting Information. The properties obtained by using this CG model were in good agreement with those obtained from AA MD simulations (see Table S2 of Supporting Information). The second step involved the development of intramolecular interactions for a PNIPAM chain and intermolecular interactions between PNIPAM and water to reproduce its LCST behavior. The chemically and structurally accurate NIPAM model was utilized to generate a PNIPAM chain with 30 monomer units (30-mer), which was obtained by placing adjacent monomers at a random angle. A 30-mer PNIPAM chain has been successfully studied to investigate the mechanism of coil-to-globule transition in many AA MD simulations by employing various FF parameters.19,40,41 The 30-mer CG PNIPAM chain was solvated with 8434 beads of 1site CG water model (see Figure 1b).35 Note, the 1-site CG model can reproduce several properties of water at various temperatures from 280 K to 370 K.35 The bonded and nonbonded interaction parameters for the backbone of PNIPAM, represented by -C2M-C2M- beads, were adopted from recently developed hydrocarbon CG model.36 The dihedral potentials were incorporated into the FF form with the goal of retaining the tacticity of PNIPAM chains (see Supporting Information Section 2 for more details).42,43 Equilibrium angle (C2M-C2M-AMP) and equilibrium dihedral angles were chosen based on the distributions obtained from mapped AA trajectory of PNIPAM simulations. During the FF parameterization, angle force constants involved with backbone atoms and dihedral force constants are set to 15 kcal/mol/radian2 and 1.0 kcal/mol, respectively. The nonbonded interactions between PNIPAM chain and water molecules were represented by 12-6 Lennard-Jones (LJ) potentials. The values of ε, the strength of interactions, between C2M-W1, AMP-W1 and ISP-W1 beads were optimized to reproduce the coil-like and globule-like conformation of PNIPAM at 290 K (below LCST) and 320 K (above LCST) (See Section 2 of Supporting Information). The value of σ, the finite distance at which the intermolecular potential between two beads becomes zero, were obtained by Lorentz-Berthelot (LB) mixing rule. Previous AA MD simulations of PNIPAM have reported, radius of gyration (Rg) value of ~16 Å and ~10 Å for PNIPAM 30-mer at ~290 K and ~320 K, respectively.19 The Rg value of ~16 Å and ~10 Å represent a coil-like and a globule-like state of PNIPAM 30-mer. In the present study, to optimize


5


Page 6 of 20

the FF parameters of PNIPAM, by using the PSO method, we set its Rgs to ~16 Å (at 290 K) and ~10 Å (320 K) as target values.19 The FF parameters were optimized by performing simulations at 290 K and 320 K temperatures simultaneously. The optimized FF parameters for PNIPAM CG model are listed in Table S3 of Supporting Information. The robustness of these FF parameters was evaluated by performing five independent simulations for 1 µs at 290 K and 320 K. The distribution of Rg obtained by averaging the last 850 ns simulation trajectory from each of these simulations (850 ns x 5 = 4250 ns) is shown in Figure 1c. This figure clearly shows a coil-like (Rg ~16 Å) and globule-like conformations (Rg ~11 Å) of 30-mer of PNIPAM below and above its LCST, respectively. Furthermore, to evaluate the ability of the new PNIPAM CG model in capturing the dependence of LCST on the PNIPAM chain length, we have generated PNIPAM chains with 5-, 18-, and 100-mers. Note, both AA MD and experiments have shown that the PNIPAM chains below 10-mers do not exhibit a coil-to-globule transition due to their short chain length.19,44–46 Two different initial configurations were generated for each chain length to perform simulations in the presence of 1-site CG water model. The simulations were carried out for 1 µs at 290 K and 320 K, which is below and above the LCST of PNIPAM, respectively. Average R g distributions obtained from both the simulations of different lengths are given in Figures 1d to 1g. Consistent with previous reported literature, PNIPAM 5-mer does not exhibit a coil-to-globule transition above its LCST.19,44,46 The peak position at both temperatures fluctuates at ~5 Å. As the chain length is increased from 5-mer to 18-mer, the PNIPAM chains exhibit a coil-like and globule-like state, below and above its LCST, respectively (Figure 1e). In the case of 30-mer and 100-mer at 290 K, there is a clear evidence of a coil-like state of PNIPAM with Rg spreading across 13 Å to 19 Å and 30 Å to 60 Å, respectively. At 320 K, the peaks in R g distributions for 30-mer and 100-mer were observed at ~11 Å and ~18 Å, respectively. Thus, the new CG PNIPAM model is able to capture the coil-to-globule transition for PNIPAM chains with different lengths.


6

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


Figure 1: (a) Mapping scheme of CG NIPAM model, (b) Schematic of CG models of PNIPAM and 1-site water,35 (c) Distribution of Rg, averaged over five independent simulations, for 30-mer at 290 K (below LCST, black lines) and 320 K (above LCST, red lines). Each CG MD simulation was performed for 1 µs in NPT ensemble and the histograms were determined for the last 850 ns (850 ns x 5 = 4250 ns). Inset shows the representative snapshots for the coil-like (black chain) and globule-like (red chain) conformations. Distribution of Rg for PNIPAM chains with (d) 5-mer, (e) 18-mer, (f) 30-mer, and (g) 100-mer at 290 K and 320 K. To evaluate the structure of water near PNIPAM chains, we calculated the radial distribution function (RDF) between 30-mer or 5-mer of PNIPAM and water beads, at two temperatures (see Supporting Information Figure S5). The PNIPAM chain with 30-mer dehydrate above the LCST which can be attributed to its globule-like state, consistent with AA results.47,48 Ahmed et. al. have predicted the dehydration time for water within 1 µs using kinetic studies on


7


Page 8 of 20

the PNIPAM collapse and our CG MD simulations of PNIPAM chains (see Figure S6 of Supporting Information) suggest comparable time-scales with these study.45 Experimentally, the coil-to-globule transition exhibited by PNIPAM is reversible i.e. when temperature is decreased from above its LCST to below the LCST, PNIPAM undergoes a globuleto-coil transition. However, Kang et al. have shown that the tightly collapsed initial configurations (a globule-like state) remain in the globule-like state below its LCST during 1 μs AA MD simulations, which suggests that these timescales could be short to observe its globule-to-coil transition.41 To further explore the robustness of our new CG model of PNIPAM in predicting such a behavior, 40 different configurations of PNIPAM with initial Rg ranging from 10 Å to 18 Å were generated. These 40 structures were binned based on their Rgs such that every 10 configurations fall in following sets: 10-12 Å (complete globule-like state) (S-I), 12-14 Å (partial globule-like state) (S-II), 14-16 Å (partial coil-like state) (S-III), and 16-18 Å (extended coil-like state) (S-IV). MD simulations of all the PNIPAM structures were carried out for 1 µs in NPT ensemble at 290 K (below LCST) and 320 K (above LCST). The Rg value for an individual conformation was determined over the last 850 ns of the total trajectory. The data obtained for individual conformations is shown in the Figure S6 of the Supporting Information. Figure 2 (a-d) shows the average distribution of Rg obtained from each set. Overall, we find that with increase in the initial Rg values, more number of PNIPAM chains showed a coil-like and a globule-like state, below and above the LCST of PNIPAM, respectively. In Figure 2 (a-d), at 290 K, the distributions are bimodal in nature since they represent the average data for 10 configurations irrespective of they were in a coil-like or a globule-like state. For S-I and S-II, the average Rg values suggest presence of PNIPAM chains in a globule-like state both below and above the LCST. In contrary, the distributions suggests that the conformations of S-III and S-IV demonstrate a clear coil-like and a globule-like state below and above the LCST of PNIPAM, respectively. These results are consistent with AA simulations reported by Kang et al.41 This further suggests that the CG model of PNIPAM, developed in this study, can accurately capture the behavior reported for AA models of PNIPAM.


8

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


Figure 2: Distribution of Rg profiles obtained by averaging the 10 different initial configurations at two temperatures: 290 K (below LCST, black lines) and 320 K (above LCST, red lines). All these ten configurations are selected based on their similar initial Rg values ranging from (a) 1012 Å (complete globule-like state) (S-I), (b) 12-14 Å (partial globule-like state) (S-II), (c) 14-16 Å (partial coil-like state) (S-III), and (d) 16-18 Å (extended coil-like state) (S-IV). Representative snapshots for each set are shown as inset to each figure.

Furthermore, to investigate the ability of the CG model in predicting the LCST of PNIPAM, we randomly selected three configurations (P1, P2, and P3) that showed an LCST transition and whose Rg was ranging from 14-18 Å (sets S-III and S-IV). If the interactions between


9


Page 10 of 20

PNIPAM and water are captured accurately then the LCST predicted by the new CG model should fall within the range observed in experiments.49 Note, it has been reported that the tacticity of PNIPAM plays an important role in determining its LCST and based on its tacticity the LCST can vary from 290 K to 305 K.49 For example, the PNIPAM chains with meso-diad content of ~45 % and ~66 % were known to show LCST at ~305 K and ~290 K, respectively.19,49 Here, the CG MD simulations were performed for 1 µs for temperatures ranging from 285 K to 320 K with an interval of 5 K. Figure 3 a-c shows the mean Rg value over last 850 ns and its standard deviation at each temperature for the PNIPAM with three different initial conformations. The LCST predicted by the three different models is ~300 K (P1), ~307 K (P2), and ~312 K (P3), which is in the comparable to the experimentally reported LCST of PNIPAM (290 K to 305 K). The differences in these LCSTs suggest that the three chains may have different tacticities. Therefore, the three simulated chains can be attributed to a meso-diad content of ~57 % (P1LCST = ~300 K), ~45 % (P2LCST = ~307 K), and less than 45 % (P3LCST = ~312 K) based on the experimental study.49 Recently, Boţan et. al. have used an autocorrelation function (ACF) of dihedral angle (ACF= , Θ being the dihedral angle at time t) to demonstrate that for an AA model of PNIPAM chain the tacticity is retained.34 Moreover, they also show that retaining such tacticity in a CG model is a very challenging task. To evaluate the ability of the newly developed PNIPAM CG model in retaining its tacticity and to further explain the differences observed in the LCST of PNIPAM, we used the ACF defined by Boţan et. al. The tacticity of the backbone dihedral angle (AMP-C2M-C2M-AMP) (See Figure 1b) for the three CG PNIPAM chains, P1, P2, and P3 at two temperatures, 290 K and 320 K was determined. As shown in Figure 3d, the dihedral angles in our CG models are highly correlated in long times, which suggests that the tacticity of PNIPAM chains is retained, similar to that of AA models and experiments.34 To the best of our knowledge, this is the first temperature-independent CG model of PNIPAM, a thermosensitive polymer, that can maintain its tacticity. This can be attributed to the presence of dihedral angles and stiffer angle force constants for C2M-C2M-C2M and C2M-C2M-AMP beads in the PNIPAM backbone. In addition, as can be seen from Figure 3d the ACF decays faster above the LCST as compared to below the LCST. This can be attributed to the coil-to-globule transition of PNIPAM chains. Similar behavior for the ACF of the dihedral angle C2M-C2M-C2M-C2M beads, that represent only the backbone, was observed (See Figure S7 of Supporting Information). We find that retaining the tacticity maintains the distance between the ISP beads present PNIPAM side-


10

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


chain. As shown in Figure 3e, the distance distribution between the alternate ISP beads shows that the ISP beads are farther in P3 conformation (whose LCST is ~312 K) than the other two conformations, namely, P1 (LCST ~300 K) and P2 (LCST ~307 K) . Decrease in the distance between ISP beads may enhance the interactions between side-chains, allowing a PNIPAM chain to undergo a coil-to-globule transition at lower temperature.

Figure 3: Evaluation of the LCST of three selected conformations from the initial Rg ranging between 14-18 Å: (a) P1LCST = ~300 K, (b) P2LCST = ~307 K, (c) P3LCST = ~312 K. (d) Autocorrelation of the dihedral angle (AMP-C2M-C2M-AMP), averaged over all the dihedrals except the edge ones. Simulation trajectories that corresponds to 290 K and 320 K are shown with solid and dashed lines, respectively. (e) Distance distribution between the alternate ISP beads for the three conformations. Color scheme used to represent the data for P1, P2, and P3 are black, red, blue, respectively.


11


Page 12 of 20

Further, we have evaluated the effect of the ability of the CG model in retaining its tacticity on its predictability of the coil-to-globule transition of PNIPAM. Specifically, we generated PNIPAM chains by placing the adjacent monomers at a certain angle ranging from 0° to 180° with an interval of 30°. The snapshots describing the different angles between the side chains are shown as insets in the Figures 4 a-g. Three initial configurations of PNIPAM chains with different backbone structures were generated to perform CG MD simulations at 290 K and 320 K for 1 μs. The 1-site CG water model was used to solvate a PNIPAM chain. The Rg for each chain was evaluated for final 850 ns of the simulation trajectory. As shown in Figures 4 a-g, a coil-to-globule transition is observed for PNIPAM chains with initial angles of 90°, 120°, and 150°. The mean values of Rg are ~18 Å and ~11 Å at 290 K and 320 K, respectively. For structures with an angle of 0°, 30°, 60°, and 180°, the PNIPAM chains were in globule-like conformation (mean Rg = ~12 Å) at both the temperatures. To gain more insights on the observed behavior, distance between the ISP beads of alternate monomers were determined for the first 150 ns from the simulation trajectories. Figure 4h illustrates the distance distribution (averaged over three different simulations) between the ISP beads of alternate monomers for the LCST (90°, 120°, and 150°) and non-LCST (0°, 30°, 60°, and 180°) PNIPAM conformations. Overall, for the conformations that showed an LCST behavior, the peak height at ~5 Å and ~11.5 Å, increases and decreases, respectively as compared to the non-LCST conformations. This may enhance the polymer-polymer interactions for non-LCST conformations and thereby PNIPAM chains form a globule-like structure below and above the LCST. Structure of water around the AMP beads of PNIPAM chains for the LCST ones (90°, 120°, and 150°) were determined at two temperatures (See Figure S8 of Supporting Information). RDF profiles suggest that the dehydration timescales for amide groups (< 1 μs) are in a reasonable agreement with the study by Ahmed et. al.45


12

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


Figure 4: Distribution of Rg for PNIPAM chains with initial angle between adjacent NIPAM monomers of (a) 0°, (b) 30°, (c) 60°, (d) 90°, (e) 120°, (f) 150°, and (g) 180° at 290 K and 320 K. All the profiles were averaged over three different initial configurations for a given angle and the corresponding values for each simulation are given in Table S4 of Supporting Information. Inset shows the one of the initial configuration. (h) Distribution of the distance between the alternate ISP beads of various PNIPAM conformations with different initial angles between the NIPAM monomers.

Next, to develop the qualitative picture of the coil-to-globule transition process, we have employed non-metric multidimensional scaling (NMDS) method, a data-driven machine-learning (ML) scheme, that reduces the dimensionality while preserving the inter-relationships of the data points (See Supplementary Information Section 5).25,50,51 A 27-dimensional vector space comprising of all the backbone dihedral angles between C2M-C2M-C2M-C2M beads of a 30-mer PNIPAM chain was utilized to understand the pathway complexity associated with the coil-toglobule transition. Note, this dihedral space was utilized because it showed differences in the ACF at two temperatures (see Figure S7 of Supporting Information). The 1 µs simulation trajectories, generated during the exploration of angle dependence between the NIPAM monomers, was considered for the NMDS analysis. Similar to other dimension reduction schemes, the axes of the


13


Page 14 of 20

NMDS data can not be simply mapped back on to the original high dimensional space. We have used symbols χ and ψ to represent the reduced dimensions from the NMDS analysis, where vertical and horizontal axes represent 2nd and 3rd dimensions respectively. To understand the mechanisms and intermediates formed during a coil-to-globule transition of PNIPAM, we have presented NMDS data in Figure 5 a, b along with the representative snapshots of PNIPAM chains. The configurations shown in Figure 5a correspond to structures with an angle of 900 between adjacent NIPAM monomers. At 320 K, Figure 5a illustrates that the PNIPAM chain goes through several metastable states, indicated by the clustering of points, before going to a globule-like state. Even below the LCST (at 290 K), PNIPAM chain exhibits several metastable coil-like states during the 1 μs simulation. Figure 5b shows the NMDS analysis of the trajectories obtained for the PNIPAM structures with different backbone configurations and an angle of 900 between adjacent NIPAM monomers. This figure suggests that a coil-to-globule transition in PNIPAM can have multiple pathways. Similar behavior has been reported for various biomolecules and hydrocarbons.25,52,53 To further understand the dependence of the LCST on the angles between adjacent monomers, we selected one trajectory of PNIPAM configuration with angle ranging from 0° to 180° for the NMDS analysis. These analyses of dihedral space suggest that the PNIPAM chains with angles 0°, 30°, 60°, and 180° transform from a coil-like to a globule-like state with 1 to 3 intermediate metastable states. On the other hand, PNIPAM conformations with 90°, 120°, and 150° angles exhibit multiple metastable states (> 3) during its coil-to-globule transition (See Supporting Information Section 5 and Figure S9). This may allow them to retain a coil-like conformation below LCST and undergo a coil-to-globule transition above the LCST.


14

Page 15 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


Figure 5: NMDS analysis showing the clustering, based on the similarity between structures in the dihedral space. (a) Evolution of a PNIPAM chain with an angle of 900 between adjacent NIPAM monomers, at two temperatures (290 K and 320 K), illustrating the LCST behavior. (b) Different pathways opted by PNIPAM chains in reaching the globule-like state at 320 K are shown for a PNIPAM chain with different backbone configurations and an angle of 900 between adjacent NIPAM monomers. Color of the points indicate the Rg value in accordance to the provided scale bar. Black and red color chains represent the snapshots that corresponds to the simulations of 290 K and 320 K, respectively. Initial PNIPAM conformation is shown in green color. Time evolution is indicated by dashed line with arrow to guide the eye.

In summary, we present a novel computational framework that integrates coarse-grained (CG) molecular dynamics (MD) simulations with a data-driven machine-learning (ML) method, which we use to study a classic thermosensitive polymer, poly(N-isopropylacrylamide) (PNIPAM). We have developed an accurate temperature-independent CG model of PNIPAM by integrating the particle swarm optimization (PSO) method with MD simulations. This is a first of its kind CG model of PNIPAM that can retain its tacticity and accurately reproduce the dependence of LCST of its tacticity. This model was thoroughly validated by predicting the LCST of PNIPAM chains with different initial conformations of the chain (different Rg as well as different angle between the adjacent monomers), varied chain length in the presence of 1-site CG water model.


15


Page 16 of 20

The analyses of 90 µs simulation trajectories clearly suggest that the new CG model could accurately reproduce the experimentally reported value of the LCST of ~305 K. Moreover, to gain insights on the pathways of the coil-to-globule transition of PNIPAM, we have utilized non-metric multidimensional scaling (NMDS) method, a data-driven ML approach, to analyze these simulation trajectories. To the best of our knowledge, this is the first study that uses a ML method to analyze the CG MD simulation trajectories of a thermosensitive polymer. This analysis reveals that the PNIPAM chains exhibiting a coil-like and a globule-like state, below and above the LCST, respectively, go through multiple metastable states. Also, a coil-to-globule transition can occur through several different pathways. This knowledge can be utilized to systematically develop a thermosensitive polymer with a desired LCST temperature and transition pathways. Moreover, the use of ML methods to unravel pathways of coil-to-globule transition demonstrates their use in studying complex systems of stimuli-sensitive polymers. This new computational framework is general and can be utilized to study complex processes like self-assembly and conformational transitions in the architectures of stimuli-sensitive polymers. ACKNOWLEDGMENT Authors would like to acknowledge Advanced Research Computing (ARC), Virginia Tech for providing the computational resources. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. S. A. D. acknowledges the Start-up Funds from Virginia Tech. Supporting Information Available: Details on the NIPAM and PNIPAM model development. Force field parameters for the developed PNIPAM CG model. Discussed the structure of water around PNIPAM chains at two temperatures. More details on NMDS analysis. REFERENCES (1) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812-7824. (2) Kenzaki, H.; Koga, N.; Hori, N.; Kanada, R.; Li, W.; Okazaki, K.-I.; Yao, X.-Q.; Takada, S. CafeMol: A Coarse-Grained Biomolecular Simulator for Simulating Proteins at Work. J. Chem. Theory Comput. 2011, 7, 1979-1989.


16

Page 17 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


(3) Rossi, G.; Monticelli, L.; Puisto, S. R.; Vattulainen, I.; Ala-Nissila, T. Coarse-Graining Polymers with the MARTINI Force-Field: Polystyrene as a Benchmark Case. Soft Matter 2011, 7, 698-708. (4) Carbone, P.; Varzaneh, H. A. K.; Chen, X.; Müller-Plathe, F. Transferability of CoarseGrained Force Fields: The Polymer Case. J. Chem. Phys. 2008, 128,064904. (5) Deshmukh, S. A.; Solomon, L. A.; Kamath, G.; Fry, H. C. Water Ordering Controls the Dynamic Equilibrium of Micelle–fibre Formation in Self-Assembly of Peptide Amphiphiles. Nature, 2016, 12367 (6) Srinivas, G.; Discher, D. E.; Klein, M. L. Self-Assembly and Properties of Diblock Copolymers by Coarse-Grain Molecular Dynamics. Nat. Mater. 2004, 3, 638-644. (7) Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable Water Model for the Coarse-Grained MARTINI Force Field. PLoS Comput. Biol. 2010, 6, e1000810. (8) Sharma, S.; Kumar, S. K.; Buldyrev, S. V.; Debenedetti, P. G.; Rossky, P. J.; Stanley, H. E. A Coarse-Grained Protein Model in a Water-like Solvent. Sci. Rep. 2013, 3, 1841. (9) Horn, J. N.; Sengillo, J. D.; Lin, D.; Romo, T. D.; Grossfield, A. Characterization of a Potent Antimicrobial Lipopeptide via Coarse-Grained Molecular Dynamics. Biochim. Biophys. Acta 2012, 1818, 212-218. (10) Hong, C.; Tieleman, D. P.; Wang, Y. Microsecond Molecular Dynamics Simulations of Lipid Mixing. Langmuir 2014, 30, 11993-12001. (11) Lee, H.; de Vries, A. H.; Marrink, S.-J.; Pastor, R. W. A Coarse-Grained Model for Polyethylene Oxide and Polyethylene Glycol: Conformation and Hydrodynamics. J. Phys. Chem. B 2009, 113, 13186-13194. (12) Winter, N. D.; Schatz, G. C. Coarse-Grained Molecular Dynamics Study of Permeability Enhancement in DPPC Bilayers by Incorporation of Lysolipid. J. Phys. Chem. B 2010, 114, 5053-5060. (13) Banerjee, P.; Roy, S.; Nair, N. Coarse-Grained Molecular Dynamics Force-Field for Polyacrylamide in Infinite Dilution Derived from Iterative Boltzmann Inversion and MARTINI Force-Field. J. Phys. Chem. B 2018, 122, 1516-1524. (14) Lu, L.; Dama, J. F.; Voth, G. A. Fitting Coarse-Grained Distribution Functions through an Iterative Force-Matching Method. J. Chem. Phys. 2013, 139, 121906. (15) Wang, Q.; Keffer, D. J.; Nicholson, D. M. A Coarse-Grained Model for Polyethylene Glycol Polymer. J. Chem. Phys. 2011, 135, 214903. (16) Prasitnok, K. A Coarse-Grained Model for Polylactide: Glass Transition Temperature and Conformational Properties. J. Polym. Res. 2016, 23, 139. (17) Peter, E. K.; Shea, J.-E.; Pivkin, I. V. Coarse kMC-Based Replica Exchange Algorithms for the Accelerated Simulation of Protein Folding in Explicit Solvent. Phys. Chem. Chem. Phys.


17


Page 18 of 20

2016, 18, 13052-13065. (18) Bejagam, K. K.; Singh, S. K.; An, Y.; Deshmukh, S. A. Machine Learnt Coarse-Grained Models. J. Phys. Chem. Lett. 2018, 9, 4667-4672. (19) Deshmukh, S. A.; Sankaranarayanan, S. K. R. S.; Suthar, K.; Mancini, D. C. Role of Solvation Dynamics and Local Ordering of Water in Inducing Conformational Transitions in poly(NIsopropylacrylamide) Oligomers through the LCST. J. Phys. Chem. B 2012, 116, 2651-2663. (20) De Silva, C. C.; Leophairatana, P.; Ohkuma, T.; Koberstein, J. T.; Kremer, K.; Mukherji, D. Sequence Transferable Coarse-Grained Model of Amphiphilic Copolymers. J. Chem. Phys. 2017, 147, 064904. (21) Chudoba, R.; Heyda, J.; Dzubiella, J. Temperature-Dependent Implicit-Solvent Model of Polyethylene Glycol in Aqueous Solution. J. Chem. Theory Comput. 2017, 13, 6317-6327. (22) Abbott, L. J.; Stevens, M. J. A Temperature-Dependent Coarse-Grained Model for the Thermoresponsive Polymer poly(N-Isopropylacrylamide). J. Chem. Phys. 2015, 143, 244901. (23) Park, W. B.; Chung, J.; Jung, J.; Sohn, K.; Singh, S. P.; Pyo, M.; Shin, N.; Sohn, K.-S. Classification of Crystal Structure Using a Convolutional Neural Network. IUCrJ 2017, 4, 486-494. (24) Ryan, K.; Lengyel, J.; Shatruk, M. Crystal Structure Prediction via Deep Learning. J. Am. Chem. Soc. 2018, 140, 10158-10168. (25) Rajan, A.; Freddolino, P. L.; Schulten, K. Going beyond Clustering in MD Trajectory Analysis: An Application to Villin Headpiece Folding. PLoS One 2010, 5, e9890. (26) Wang, J.; Gayatri, M.; Ferguson, A. L. Coarse-Grained Molecular Simulation and Nonlinear Manifold Learning of Archipelago Asphaltene Aggregation and Folding. J. Phys. Chem. B 2018, 122, 6627-6647. (27) Izenman, A. J. Introduction to Manifold Learning. WIREs Comp Stat 2012, 4, 439–446. (28) Jain, K.; Vedarajan, R.; Watanabe, M.; Ishikiriyama, M.; Matsumi, N. Tunable LCST Behavior of poly(N-Isopropylacrylamide/ionic Liquid) Copolymers. Polym. Chem. 2015, 6, 6819–6825. (29) Guan, Y.; Zhang, Y. PNIPAM Microgels for Biomedical Applications: From Dispersed Particles to 3D Assemblies. Soft Matter 2011, 7, 6375–6384. (30) Gandhi, A.; Paul, A.; Sen, S. O.; Sen, K. K. Studies on Thermoresponsive Polymers: Phase Behaviour, Drug Delivery and Biomedical Applications. Asian J. Pharm. Sci. 2015, 10, 99– 107. (31) Zhernenkov, M.; Ashkar, R.; Feng, H.; Akintewe, O. O.; Gallant, N. D.; Toomey, R.; Ankner, J. F.; Pynn, R. Thermoresponsive PNIPAM Coatings on Nanostructured Gratings for Cell Alignment and Release. ACS Appl. Mater. Interfaces 2015, 7, 11857–11862.


18

Page 19 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


(32) Wang, Z.; He, X. Developing Coarse-Grained Force Fields for PNIPAM Single Chain from the Atomistic Model. In 2009 3rd International Conference on Bioinformatics and Biomedical Engineering; 2009; 1–4. (33) Abbott, L. J.; Stevens, M. J. A Temperature-Dependent Coarse-Grained Model for the Thermoresponsive Polymer poly(N-Isopropylacrylamide). J. Chem. Phys. 2015, 143, 244901. (34) Boţan, V.; Ustach, V. D.; Leonhard, K.; Faller, R. Development and Application of a CoarseGrained Model for PNIPAM by Iterative Boltzmann Inversion and Its Combination with Lattice Boltzmann Hydrodynamics. J. Phys. Chem. B 2017, 121, 10394–10406. (35) Bejagam, K. K.; Singh, S.; An, Y.; Berry, C.; Deshmukh, S. A. PSO-Assisted Development of New Transferable Coarse-Grained Water Models. J. Phys. Chem. B 2018, 122, 1958-1971. (36) An, Y.; Bejagam, K. K.; Deshmukh, S. A. Development of New Transferable Coarse-Grained Models of Hydrocarbons. J. Phys. Chem. B 2018, 122, 7143–7153. (37) Bejagam, K. K.; Singh, S.; Deshmukh, S. A. Development of Non-Bonded Interaction Parameters between Graphene and Water Using Particle Swarm Optimization. J. Comput. Chem. 2017, 39, 721-734. (38) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671–690. (39) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. (40) Chiessi, E.; Paradossi, G. Influence of Tacticity on Hydrophobicity of Poly(NIsopropylacrylamide): A Single Chain Molecular Dynamics Simulation Study. J. Phys. Chem. B 2016, 120, 3765–3776. (41) Kang, Y.; Joo, H.; Kim, J. S. Collapse–Swelling Transitions of a Thermoresponsive, Single Poly(N-Isopropylacrylamide) Chain in Water. J. Phys. Chem. B 2016, 120, 13184–13192. (42) Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A.; Kremer, K. Hierarchical Modeling of Polystyrene: From Atomistic to Coarse-Grained Simulations. Macromolecules 2006, 39, 6708–6719. (43) Li, C.; Shen, J.; Peter, C.; van der Vegt, N. F. A. A Chemically Accurate Implicit-Solvent Coarse-Grained Model for Polystyrenesulfonate Solutions. Macromolecules 2012, 45, 2551– 2561. (44) Abbott, L. J.; Tucker, A. K.; Stevens, M. J. Single Chain Structure of a poly(NIsopropylacrylamide) Surfactant in Water. J. Phys. Chem. B 2015, 119, 3837–3845. (45) Ahmed, Z.; Gooding, E. A.; Pimenov, K. V.; Wang, L.; Asher, S. A. UV Resonance Raman


19


Page 20 of 20

Determination of Molecular Mechanism of poly(N-Isopropylacrylamide) Volume Phase Transition. J. Phys. Chem. B 2009, 113, 4248–4256. (46) Deshmukh, S. A.; Sankaranarayanan, S. K. R. S.; Mancini, D. C. Vibrational Spectra of Proximal Water in a Thermo-Sensitive Polymer Undergoing Conformational Transition across the Lower Critical Solution Temperature. J. Phys. Chem. B 2012, 116, 5501–5515. (47) Deshmukh, S. A.; Li, Z.; Kamath, G.; Suthar, K. J.; Sankaranarayanan, S. K. R. S.; Mancini, D. C. Atomistic Insights into Solvation Dynamics and Conformational Transformation in Thermo-Sensitive and Non-Thermo-Sensitive Oligomers. Polymer 2013, 54, 210–222. (48) Futscher, M. H.; Philipp, M.; Müller-Buschbaum, P.; Schulte, A. The Role of Backbone Hydration of Poly(N-Isopropyl Acrylamide) Across the Volume Phase Transition Compared to Its Monomer. Sci. Rep. 2017, 7, 17012. (49) Ray, B.; Okamoto, Y.; Kamigaito, M.; Sawamoto, M.; Seno, K.-I.; Kanaoka, S.; Aoshima, S. Effect of Tacticity of Poly (N-Isopropylacrylamide) on the Phase Separation Temperature of Its Aqueous Solutions. Polym. J. 2005, 37, 234. (50) Taguchi, Y.-H.; Oono, Y. Nonmetric Multidimensional Scaling As a Data-Mining Tool: New Algorithm and New Targets: Applications to Chemical Reaction Dynamics in Complex Systems. In Geometric Structures of Phase Space in Multidimensional Chaos; Toda, M., Komatsuzaki, T., Konishi, T., Berry, R. S., Rice, S. A., Eds.; Advances in Chemical Physics; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2005, 315–351. (51) Pereda, M.; Estrada, E. Machine Learning Analysis of Complex Networks in Hyperspherical Space. arXiv [physics.soc-ph], 2018. (52) Ferguson, A. L.; Panagiotopoulos, A. Z.; Kevrekidis, I. G.; Debenedetti, P. G. Nonlinear Dimensionality Reduction in Molecular Simulation: The Diffusion Map Approach. Chem. Phys. Lett. 2011, 509, 1–11. (53) Ferguson, A. L.; Panagiotopoulos, A. Z.; Debenedetti, P. G.; Kevrekidis, I. G. Systematic Determination of Order Parameters for Chain Dynamics Using Diffusion Maps. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 13597-13602.


20

Machine-Learning Enabled New Insights into the Coil-to-Globule

Recommend Documents