Sequence transferable coarse-grained model of amphiphilic copolymer

TL;DR

This paper develops a transferable coarse-grained model for amphiphilic copolymers that captures their multi-scale behavior and is validated against experimental data, enabling better simulation of complex polymer systems.

Contribution

A novel segment-based coarse-grained model for amphiphilic copolymers that is transferable across different sequences and consistent with experimental observations.

Findings

Model accurately reproduces global conformational behavior.

Transferable across various copolymer sequences.

Aligns well with experimental data.

Abstract

Polymer properties are inherently multi-scale in nature, where delicate local interaction details play a key role in describing their global conformational behavior. In this context, deriving coarse-grained (CG) multi-scale models for polymeric liquids is a non-trivial task. Further complexities arise when dealing with copolymer systems with varying microscopic sequences, especially when they are of an amphiphilic nature. In this work, we derive a segment-based generic CG model for amphiphilic copolymers consisting of repeat units of hydrophobic (methylene) and hydrophilic (ethylene oxide) monomers. The system is a simulation analogue of polyacetal copolymers [Samanta et al., Macromolecules 49, 1858 (2016)]. The CG model is found to be transferable over a wide range of copolymer sequences and also to be consistent with existing experimental data.

Tables

Table 1. Table 1: Copolymer conformations with different sequences (as shown in the first four columns). In the fifth column we present methylene mole fraction x m subscript 𝑥 𝑚 x_{m} along the copolymer backbone, which is defined as x m = ( n 1 + n 2 ) / ( n 1 + m 1 + n 2 + m 2 ) subscript 𝑥 𝑚 subscript 𝑛 1 subscript 𝑛 2 subscript 𝑛 1 subscript 𝑚 1 subscript 𝑛 2 subscript 𝑚 2 x_{m}=(n_{1}+n_{2})/(n_{1}+m_{1}+n_{2}+m_{2}) . The data is shown for T = 320 𝑇 320 T=320 K. The transition temperature in our reference all-atom simulations is about 20 K lower than the experimental transition temperature. Therefore, we present the extrapolated polymer conformations from experiments, with the transition temperature shifted lower by 20 K. In the last column, we present corresponding polymer conformation from CG simulations. Note that the polymer conformations are identified by calculating their single chain form factor S ​ ( q ) 𝑆 𝑞 S(q) , which shows S ​ ( q ) ∼ q − 5 / 3 similar-to 𝑆 𝑞 superscript 𝑞 5 3 S(q)\sim q^{-5/3} for an expanded chain and S ​ ( q ) ∼ q − 4 similar-to 𝑆 𝑞 superscript 𝑞 4 S(q)\sim q^{-4} for a collapsed globule.

$n_1$	$m_1$	$n_2$	$m_2$	$x_m$	Experiment [6]	Simulation
4	0	2	0	1.00	Turbid	Globule
4	0	2	1	0.86	Turbid	Globule
4	0	2	2	0.75	Turbid	Globule
4	0	2	3	0.67	Turbid	Globule
4	0	2	4	0.60	Turbid	Globule
4	0	2	5	0.54	Turbid	Globule
2	1	2	1	0.67	Turbid	Globule
2	1	2	2	0.57	Turbid	Globule
2	1	2	3	0.50	Clear	Lamellar
2	1	2	4	0.44	Clear	Expanded
2	1	2	5	0.40	Clear	Expanded
2	2	2	1	0.57	Clear	Globule
2	2	2	2	0.50	Clear	Expanded
2	2	2	3	0.44	Clear	Expanded
2	2	2	4	0.40	Clear	Expanded
2	2	2	5	0.36	Clear	Expanded

Full text

I Introduction

Polymer processing requires a detailed understanding of their structure-property relationships, which not only is fundamentally challenging but also has a wide range of applications ranging from physics to biology [1, 2]. In this context, smart responsive polymers serve as excellent candidates. A polymer is commonly known as “smart responsive” when a small change in external stimulus can drastically change its structure, function and/or stability. Furthermore, in these systems, where the relevant energy scale is of the order of the thermal energy $k_{\rm B}T$, large conformational and local compositional fluctuations play a delicate role dictating their properties [3, 4]. When the external stimulus is temperature $T$, these polymers are referred to as thermoresponsive polymers [5, 6]. One of the interesting classes of thermoresponsive polymers is those that remain expanded at low $T$ and collapse into compact objects when $T>T_{\ell}$, where $T_{\ell}$ is referred to as a lower critical solution temperature (LCST). The properties of LCST polymers are usually dictated by hydrogen bonding of water molecules with polymer segments, which break down when $T>T_{\ell}$. During this process, when hydrogen hydrogen bonds are broken, water molecules are expelled from near the polymer structure and thus gain a large amount of translational entropy. Therefore, simply speaking the polymer chain will collapse whenever the total translational entropy gained by water molecules is larger than the conformational entropy lost by the polymer upon collapse [7].

For a given homopolymer structure there exists a characteristic $T_{\ell}$ [8, 9, 10]. This $T_{\ell}$ can, however, be tuned by introducing more hydrophobic or more hydrophilic units along the homopolymer backbone [6, 11, 12, 13, 14, 15]. In particular, $T_{\ell}$ increases (decreases) with increasing hydrophilic (hydrophobic) units. Moreover, the changes in $T_{\ell}$ are usually difficult to predict and non-linear with changing copolymer sequences [11, 12, 13]. However, in a recent work, based on the acetal linkage of methylene and ethylene oxide units, it has been shown that $T_{\ell}$ can be linearly tuned by adjusting the fractions of methylene and ethylene oxide monomers along the polymer backbone [6]. In addition to predictability with changes in sequence, the acetal-based copolymers are advantageous due to their biodegradable properties. Moreover, these polymers show strong chain length effects dependent on end-group nature (see Fig. 5 in Ref. [6]). For example, it has been shown that a molecular weight $M_{\rm w}\sim 10^{4}$ g/mol is necessary to avoid chain length effects, which corresponds to an end-to-end distance $R_{\rm ee}\sim 10$ nm for a poly(ethylene oxide). Within a simulation setup this would require a box size of at least 15 to 20 nm. For a water box, with number density of 32 water molecules per cubic nm, more than $10^{5}$ water molecules would be needed (or equivalent of $\sim 3\times 10^{5}$ atoms). This poses a serious challenge for all-atom simulations of these systems and motivates the development of accurate, lower resolution CG model for further theoretical investigations. While there are studies in deriving implicit and explicit solvent CG models of PEO [16, 17, 18], in this work we derive a sequence transferable, systematic CG model of an analogue of polyacetal that represents each copolymer units with single CG beads. We employ the model to study different copolymer architectures that are complementary to the earlier experimental results [6]. Additionally, we also propose a broad range of molecular morphologies, going beyond available experimental data. Note that here we discuss a case where hydrophobicity (or hydrophilicity) is tuned along the backbone and not added as a side group. However, there are cases, such as poly(N-isopropylmethacrylamide) where an extra methylene group added as a side group to poly(N-isopropylacrylamide) increases $T_{\ell}$ [19]. The possible reason behind an added hydrophobic group increases $T_{\ell}$ can be attributed to the solvent structure around the chain [15].

II Model and method

We start by commenting on the underlying reference all-atom simulations. For this purpose, the chemical structures are taken to be the same as in experiment [6]. In Fig. 1(a) we show a typical chemical configuration of the polyacetal copolymer architecture.

In this amphiphilic structure, $T_{\ell}$ is tuned by changing the numbers of methylene ($n_{1}$ and $n_{2}$) and ethylene oxide ($m_{1}$ and $m_{2}$) units. For computational simplicity of the computing, we have chosen the repeat unit of the chain as $N=5$ (see Fig. 1(b)). Here we consider a copolymer structure represented with $n_{1}=4$, $m_{1}=0$, $n_{2}=2$, and $m_{2}=3$. Note that the molecular weight of the simulated chain presented in Fig. 1(b) is almost one order of magnitude smaller than the $M_{\rm w}$ needed to avoid strong chain length effects [6]. However, it is noteworthy that the system size effects in the experiments are associated with the end group effects. Therefore, to avoid the system size effect and to make a reasonable estimate of $T_{\ell}$, we have terminated the ends of copolymer with an inert CH3 groups, as shown in Fig. 1(b). As will be discussed in the latter part of this manuscript, this approximation allows us to obtain a reasonable $T_{\ell}$, while not attempting to make any claims on the system size effects in simulations. Moreover, this small chain length can not be used for any structural predictions that only happen for the longer chains.

All-atom molecular dynamics (MD) simulations are performed using the GROMACS package [20]. The temperature is varied from 290 K to 320 K using velocity rescaling with a coupling constant 0.5 ps [21]. Each of these simulations are performed for 50 ns production runs, which is at least one order of magnitude larger than the longest relaxation time. The average is taken over the last 10 ns of MD data. The electrostatics are treated using Particle Mesh Ewald [22]. The interaction cut-off for non-bonded interactions is chosen as 1.0 nm. The simulations are performed with a constant pressure ensemble, where the the pressure is controlled using a Berendsen barostat [23] with a coupling time of 0.5 ps and 1 atm pressure. The time step for the simulations is set to 2 fs and the equations of motion are integrated using the leap-frog algorithm. LINCS algorithm is used to constraint all bond vibrations [24].

We find the gyration radius $R_{\rm g}$ to be $1.33\pm 0.15$ nm (for $T=290$ K), $1.38\pm 0.16$ nm (for $T=300$ K), $0.97\pm 0.04$ nm (for $T=310$ K), and $0.98\pm 0.09$ nm (for $T=320$ K) for the system presented in Fig. 1(b). There is a reasonably sharp change in $R_{\rm g}$ between 300K and 310 K, suggesting $T_{\ell}$ to be between 300 K and 310 K. The experimental phase diagram, however, suggests $T_{\ell}\sim 320$ K for the same copolymer sequence. Therefore, our all-atom data underestimates $T_{\ell}$ by about 10$-$20 K.

Non-bonded CG potentials are derived using the structure-based techniques for solutions [25, 26, 27]. To derive the CG model, which is sequence transferable, we make two assumptions:

1. Upon changing composition and sequencing along the backbone, there are no cross-correlations between different monomer units. In this context, experiments have shown that $T_{\ell}$ changes linearly with changing fractions of hydrophobic or hydrophilic units [6]. Therefore, presenting a situation where a segment based coarse-graining (only building CG models based on separate simulations of different monomer units) can be performed, which can later be incorporated into a polymer chain with different sequences.

In this context, a generic study of copolymer sequences, within a mean-field picture, have shown linear or non-linear pair interpolation depending of the interaction parameters [14]. 2) The acetal linker (represented by red in Fig. 1) only contributes to a negligible shift in $T_{\ell}$. Therefore, we do not incorporate an acetal linker in our simulations. We note that one might use a hard sphere type representation of the acetal unit, although this is beyond the motivation of the present study. The step-by-step procedure for building the CG model is: For the non-bonded interaction, we perform all-atom simulation of different monomer unitss solvated in water, with the mapping scheme presented in Fig. 2(a).

III Results and discussions

The resultant interactions potentials [28], shown in Fig. 3, are used to simulate the whole range of copolymer sequences. Note that the derived potentials from the monomer level is used in simulations of polymer chain with one-four exclusion of bonded neighbors. Importantly, we derive only one set of $u(r)$ that we will use for a wide range of copolymer sequences. Bonded interactions are obtained by Boltzmann inversion of bonded, angle and dihedral distributions. CG simulations are also performed in GROMACS package [20], where the interaction cut-off is chosen as 1.8 nm and the time step is chosen as 2 fs. Simulations are performed for a 100 ns long CG MD trajectory.

We start by investigating the conformation of alkane and poly(ethylene oxide) chains in water. For this purpose, the molecular weight is chosen as $M_{\rm w}=10^{4}$ g/mol. The CG model reproduces the expected conformations for both these chains, i.e., an expanded chain for poly(ethylene oxide) with $R_{\rm g}=3.91\pm 0.46$ nm, which is consistent with the previous simulations [29] and a collapsed structure for the alkane chain with $R_{\rm g}=1.18\pm 0.05$ nm. The structures are characterized by their single chain structure factor $S(q)$ shown in Fig. 4. For example, expanded chain shows a scaling law $q^{-1/\nu}$ with $\nu=3/5$ being the Flory’s exponent for a good solvent chain (see black curve in Fig. 4). When a polymer collapses, its conformation is well described by a hard sphere scattering function with a scaling law $q^{-4}$ (see red curve in Fig. 4) [7].

Now we focus on the main theme of this work, i.e., to study the conformational behavior of complex amphiphilic structures and their comparison to known experimental data obtained from the polyacetal system [6]. In Table 1 we summarize data from experimental copolymers and its comparison to these CG simulations. Note that turbid solutions found by experiment indicate collapsed structures, while clear solutions are associated with expanded conformations. Comparing experimental and simulation data in the last two columns, we find that the polymer conformation is reasonably well reproduced by the CG model, except for two cases. We would also like to point out that for methylene mole fraction $x_{m}>50\%$, chains are always collapsed as shown in the column five of Table 1.

We also would like to draw attention to the fact that it is generally challenging to incorporate composition (or sequence) and temperature transferability in a CG model. Table 1 demonstrates that our CG model is well transferable over a wide range of copolymer sequences. Further testing reveals that the CG model only reproduces reliable structures for this particular thermodynamic state of parameterization, i.e., 320 K. This is not surprising given that the many-body potential of mean force (PMF) is state point-dependent. However, there are examples where the many-body features are weak and thus the CG model can be $T$ transferable. Moreover, hydrogen bonding nature of the interactions, as in the case of polyacetal, adds up an additional complexity, making the many-body effect more relevant and thus CG potential become less transferable.

The advantage of a segment-transferable CG model, which originate because of the lack in cross-correlation between different monomer units over large length scales along the polymer backbone, is not only that it allows a rather consistent comparison with the experiments [6], but that is also enables structural predictions to be made for several more macromolecular architectures. Therefore, the CG model presented here can be viewed as a molecular toolbox to investigate the properties of different polymer architectures for advanced functional uses [28]. In this context, it has been previously predicted that amphiphilic copolymers can exhibit interesting structures [30, 31], which was later studied by generic simulations [32] and Monte Carlo simulations [33]. These complex structures are highly interesting for biomedical applications, such as drug delivery and materials for tissue engineering [1, 33]. Therefore, to investigate a broad range of conformational properties of polyacetal-based systems, we have constructed a set of 192 copolymer configurations in water with varying amphiphilic sequences. Each of these simulations were performed for 100 ns in CG units, generating a total of about 15 $\mu$s of CG MD data.

In Fig. 5 we present three sets of representative phase diagrams for the copolymer representing with varying $n_{1}$, $n_{2}$, $m_{1}$, and $m_{2}$. It can be seen that depending on the sequences, we observe a variety of copolymer configurations. For example, we find flower, micellar, multiple-flower, lamellar, and/or worm-like micellar structures. This presents a fully flexible and versatile molecular tool box for the simulation of these amphiphilic systems. In addition to the 192 configurations presented here many additional configurations may be generated by employing distinct combinations of hydrophobic and hydrophilic units. Furthermore, the applicability of this CG model is not only restricted to the linear amphiphilic chains, but may also be useful to study the solvation and/or aggregation of branched or brush like copolymers [34].

IV Conclusions

We have presented a protocol to obtain a sequence-transferable coarse-grained (CG) model for amphiphilic copolymer architectures. Our CG model was derived from a segment (monomer) based level and then was translated into longer chain simulations with different amphiphilic sequences. The derived CG model was validated by one-to-one comparison with experimental data [6]. Additionally, we also make several structural predictions that go beyond the polymers synthesized in experiments. While the derived CG model is sequence transferable, it shows poor temperature transferability. This is due to the hydrogen bonded nature of the underlying interaction details, which leads to a complex many-body effects that can not be captured within structure-based CG model. Here, however, parameterization of a CG model at different temperatures $T$ can be performed to obtain a $T$ dependent interaction capturing correct nature of hydrogen bonded interaction [35, 37]. Furthermore, the segment transferability is important because one set of CG potentials can describe, predict, and validate many amphiphilic polymer architectures. Therefore, the potentials presented here can be viewed as a molecular toolbox for a wide range of alkane and ethylene oxide based architectures.

V Acknowledgement

We thank Burkhard Dünweg and Carlos M. Marques for countless interesting discussions and Tiago E. de Oliveira for the help with the $\mathcal{C}-$IBI scripts implemented in the VOTCA package [36] and all atom force field. We further acknowledge Joseph Rudzinski, Nancy Carolina Forero-Martinez, and Torsten Stühn for critical reading of the manuscript. C.D.S., P.L. and J.T.K. wish to acknowledge support by the National Science Foundation under grant DMR-1206191. D.M. thank Joachim Dzubiella for discussions related to his manuscript draft on CG model of poly(ethyle oxide) homopolymer in water [37].