PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Polymer (Guildf). Author manuscript; available in PMC 2010 July 31.
Published in final edited form as:
Polymer (Guildf). 2009 July 31; 50(16): 4139–4149.
doi:  10.1016/j.polymer.2009.06.055
PMCID: PMC2742384
NIHMSID: NIHMS128961

A systematic procedure to build a relaxed dense-phase atomistic representation of a complex amorphous polymer using a coarse-grained modeling approach

Abstract

A systematic procedure has been developed to construct a relaxed dense-phase atomistic structure of a complex amorphous polymer. The numerical procedure consists of (1) coarse graining the atomistic model of the polymer into a mesoscopic model based on an iterative algorithm for potential inversion from distribution functions of the atomistic model, (2) relaxation of the coarse grained chain using a molecular dynamics scheme, and (3) recovery of the atomistic structure by reverse mapping based on the superposition of atomistic counterparts on the corresponding coarse grained coordinates. These methods are demonstrated by their application to construct a relaxed, dense-phase model of poly(DTB succinate), which is an amorphous tyrosine-derived biodegradable polymer that is being developed for biomedical applications. Both static and dynamic properties from the coarse-grained and atomistic simulations are analyzed and compared. The coarse-grained model, which contains the essential features of the DTB succinate structure, successfully described both local and global structural properties of the atomistic chain. The effective speedup compared to the corresponding atomistic simulation is substantially above 102, thus enabling simulation times to reach well into the characteristic experimental regime. The computational approach for reversibly bridging between coarse-grained and atomistic models provides an efficient method to produce relaxed dense-phase all-atom molecular models of complex amorphous polymers that can subsequently be used to study and predict the atomistic-level behavior of the polymer under different environmental conditions in order to optimally design polymers for targeted applications.

1. INTRODUCTION

The family of biodegradable polymers referred to as tyrosine-derived polyarylates (TDPAs) are currently under investigation for possible use in a wide range of medical devices such as vascular grafts, drug delivery vehicles, and tissue-engineering scaffolds [1-4]. Candidate members of this polymer family are being synthesized by the Kohn laboratory at Rutgers University through the use of combinatorial chemistry techniques to prepare a series of structurally related diblock TDPAs derived from monomers consisting of a tyrosine-derived diphenol and a diacid, each with various pendant groups [5-6]. The basic structure of the diphenol consists of a unit of “desaminotyrosine” and a unit of L-tyrosine alkyl ester, linked together via a peptide bond. These polymers are devised to provide desirable engineering properties with good processibility.

As part of the development of these biodegradable TDPAs, which break down by hydrolysis, it is important to understand how water diffuses into the material and partitions itself around the functional groups of the polymer as a function of its 3-D structure. Protein adsorption behavior to the surface of these types of polymers is also of interest because the bioactive states of proteins that adsorb onto a polymer surface following implantation in the body govern the cellular response to the material [7]. While combinational synthesis and testing combined with neural network analysis approaches have been applied to TDPAs to identify correlations between the structure of these polymers and their interactions with water and proteins [8-12], these types of analyses are limited in that they do not provide direct insights into the cause-and-effect relationships at the atomistic level that determine the behavior of the system. To investigate system behavior at this level, equilibrated, dense-phase molecular models of both bulk and surface structures of the polymer are required that can be used to simulate the diffusion and partitioning of water in the polymer as well as the adsorption behavior of proteins to the surface. However, the development of equilibrated, dense-phase molecular models of TDPA polymers is problematic because of their complex, branched monomeric structures combined with the relatively high molecular weight of the polymer, which is typically around 100 kDa. Given the fact that the longest relaxation time of an entangled polymer network of chain length N scales at least as N3, leading to more than N4 in CPU time for adequate system sampling [13], the use of standard all-atom molecular simulation methods to equilibrate such a molecular system is prohibitively expensive.

To circumvent this problem, coarse-grained (CG) methods combined with forward and reverse mapping techniques between all-atom and CG molecular models can be used as an approach to develop equilibrated, dense-phase models of TPDA polymers. These methods provide a means of greatly reducing the number of degrees of freedom (DOF) in the system while keeping the description of macroscopic properties at the necessary level of accuracy and detail. To map an atomistic model onto a CG model, several atoms are grouped together into relatively simple “super-atoms” and the interaction energies between two super-atoms are obtained by the application of an optimization procedure that reproduces the structural distributions obtained from atomistic simulations. The super-atom in a CG model can be either spherical, as that used in most of the published papers [14-22], or ellipsoidal by the use of anisotropic Gaussian potentials [23]. CG simulation methods have been successfully applied in various simpler macromolecular systems, such as the melts of polyethylene and cis-polybutandiene [16], various polycarbonate melts [17-19], diblock copolymer poly(styrene-b-butadiene) [20], and diblock copolymer assemblies [21, 22]. With the CG methods, the bulk structure of polymers or their assemblies can be equilibrated in reasonable time frames, and the atomistic bulk structures can be rebuilt from superpositions of atomistic counterparts on the corresponding coarse grained entities [15, 18].

Based on these successful developments and applications of CG models to address polymeric systems, we have developed and applied a novel, efficient, systematic approach to coarse grain the atomistic structure of a complex amorphous long-chain polymer, relax the structure of polymer chains based on the CG model, and the rebuild the atomistic structure from the relaxed CG structure to create dense-phase, all-atom molecular models of the amorphous polymer. In this paper, we demonstrate the implementation of these methods for poly(DTB succinate), which is a member in the TDPA family of polymers.

2. METHODOLOGY

The methods that we used for the development of the relaxed, dense-phase, all-atom model of poly(DTB succinate) are presented below in three subsections. In Section 2.1, we present the equations and approach that we used for developing the CG parameters for poly(DTB succinate). In Section 2.2, we present the methods that we used to conduct hybrid molecular dynamics (MD) and Monte Carlo (MC) simulations of all-atom segments of the polymer chain from which CG parameters were determined by the application of the equations presented in Section 2.1. Then, in Section 2.3, we present our methods for the relaxation and reverse-mapping steps that are applied for final model development.

2.1. Coarse - Grain Parameterization of Poly(DTB succinate)

Poly(DTB succinate) is characterized by a pendant butyl chain in the diphenol component of the diblock copolymer (Fig. 1A). In this work, a 7:1 mapping scheme was developed to coarse grain the DTB succinate polymer (Fig. 1B). In this mapping scheme, each monomeric repeat unit of the polymer was represented by seven spherical beads, which correspond to the two phenyl-formate (PhF), the two ethyl (Et1 and Et2), the n-ethyl-formamide (nEF), the methyl-formate (MeF), and the propyl (Pro) subunits. The basic criteria that we followed to decide which groups of atoms to combine into individual beads were: (1) The structural data of the small molecular segments represented by each bead had to be readily available for similar small molecule analogs so they could be used to establish the van der Waals (vdW) interactions for the CG segments (see Table 2), (2) the molecular segments had to be net neutral in their electrostatic charge, (3) the molecular segment effectively had to represent the local rigidity of the polymer chain (i.e., don’t combine stiff and compliant polymer segments into one bead), and (4) the molecular segments had to be able to effectively represent differences in the local chain structure. For instance, the two different types of ethyl subunits were represented because the chemical environment around Et1 and Et2 is different, leading to different distributions of structural factors. Once the bead segments were decided upon, a mapping point was then defined for each bead at its center of mass. The structural factors of the CG model of DTB succinate are summarized in Table 1.

Fig. 1Fig. 1
(A). Schematic representation of the CG model for DTB succinate. Each chemical repeat unit is represented by seven beads, corresponding to the n-ethyl-formamide (nEF), the methyl-formate (MeF), the propyl (Pro), the two phenyl-formate (PhF) and the two ...
Table 1
Summary of the atom connections in the CG model of DTB succinate.
Table 2
Liquid state properties of the compounds that have similar structures to the super-atoms in the CG model of DTB succinate.

The corresponding CG force field is developed using the CHARMM molecular simulation package [24, 25] with PCFF force field parameters [26-29], in which the all-atom model total force field energy (E) is calculated by

E=Ebonded+Enon-bonded,
(1)

where Ebonded and Enon-bonded are bonded and nonbonded potential terms, respectively, which are decomposed as,

Ebonded=Estretch+Ebend+Etorsion+EcrossEnon-bonded=EvdW+ECoulomb,
(2)

where Estretch, Ebend and Etorsion are the bond-length-stretching, angle-bending, and bond-rotation energies, respectively; the cross term, Ecross, accounts for the coupling between individual bonded interaction; EvdW accounts for the excluded volume repulsive as well as the intermolecular attractive forces between non-bonded atoms; and ECoulomb is the electrostatic potential. The CG force field is developed in a similar manner with a few modifications. Specifically, the Ecross term is not considered and the electrostatic interactions between non-bonded super-atoms are reasonably ignored since all CG beads are neutral. Therefore, the potential terms that need to be parameterized include Estretch, Ebend, Etorsion, and EvdW. For each term, the parameterization was conducted based on the iterative Boltzmann inversion of the corresponding atomistic distribution functions [30]. During the iterative procedure for a particular distribution, the rest of the potentials were kept constant. In dense systems, individual distributions usually depend on the full set of potentials through higher-order correlations, requiring that individual potentials must be readjusted after other interactions are changed. In practice, the iteration is usually started with the potentials that are least affected by changes of others. According to the suggestion of Reith et al. [31], the potential terms were therefore optimized in the order of their relative strength, with parameterization systematically developed in the order of EstretchEbendEvdWEtorsion. Global re-adjustment of the parameters for each bonded term was conducted after the vdW parameters were determined.

The Boltzmann-inversion method performs potential inversion from a set of known distributions of structural parameters (e.g., from the atomistic simulation) to extract effective CG potentials. For our system, the potentials calculated from the Boltzmann-inversion methods must reproduce the distributions of structural parameters including six different bonds, eight angles, nine dihedral angles, and five intermolecular radial distribution functions extracted from reference atomistic simulations. For poly(DTB succinate ), the four types of potentials are calculated as follows.

2.1.1. Bond Stretching

Following the work of Tschöp et al [32], the normalized distribution of bond length is calculated by

P(l)=Zlp(l)l2=exp[EstretchkBT],
(3)

where l is the distance between successive beads, Zl is a normalizing factor, l2 is a weighting factor to account for spherical volume element, and P(l) and p(l) are normalized and unnormalized distribution functions of l, respectively. From the analysis of the fluctuation of the distance between the center of mass of successive super-atoms along the chain, we found that the normalized distributions of bond length can be effectively represented by a double Gaussian function, which is expressed as

P(l)=A1w1exp(2(llc1)2w12)+A2w2exp(2(llc2)2w22)=exp[EstretchkBT],
(4)

where A1, A2, lc1, lc2, w1, and w2 are the parameters obtained through fitting, kB is the Boltzmann factor, and T is absolute temperature. Then, taking the logarithm of both sides of the equation and performing the Boltzmann-inversion, one obtains

Estretch=kBT{Kb1(llc1)2ln[1+Kb3exp[Kb1(llc1)2Kb2(llc2)2]]},
(5)

where the energetic parameters are defined as: Kb1=2/w12, Kb2=2/w22, and Kb3=w1 A2/w2A1. The parameter Kb3 accounts for the relative position of the two Gaussian peaks.

2.1.2. Bond Angle Bending

Following the treatment of Tschöp et al [32], the distributions of bond angle were weighted by a factor sin(θ) and renormalized by a factor Zθ. The normalized distribution was calculated by

P(θ)=Znp(θ)sin(θ)=exp[EbendkBT],
(6)

where θ is the angle between successive bonds, and P(θ) and p(θ) are normalized and unnormalized distribution functions of θ, respectively. The distributions of bond angles between successive CG bonds, P(θ), were also fitted with double Gaussian functions, and then were Boltzmann-inverted to calculate Ebend; i.e.,

Ebend=kBT{Kθ1(θθc1)2ln[1+Kθ3exp[K1(θθc1)2Kθ2(θθc2)2]]},
(7)

where Kθ1, Kθ2 and Kθ3 are fitted parameters determining the bending strength. They are defined similarly as the parameters in Eqn. (5).

2.1.3. Dihedral Angle

In PCFF, the torsional potential is calculated by the Fourier progression formula, i.e.,

Etorsion(ϕ)=n=13Kn[1+cos(nϕδn)],
(8)

where ϕ is the dihedral angle, Ki and δi (i = 1, 3) are force constants and phase angles, respectively. The torsional potential was calculated similarly in the CG model developed in this work. Based on the target distributions of the nine types of dihedral angles calculated from atomistic model, the parameters Ki and δi in the CG model were optimized so that the target probability distributions of dihedral angle could be reproduced by CG simulation.

2.1.4. Nonbonded Interactions

A novel scheme was developed and implemented in this work to parameterize nonbonded interactions. Only the vdW term is optimized in our current CG approach. In the atomistic model of the polymer, the PCFF force field represents vdW interactions between a pair of non-bonded atoms, i and j, by using a Lennard-Jones 6 - 9 potential; i.e.,

Evdw=εij[2(σijrij)93(σijrij)6],
(9)

where εij is the depth of the potential well, σij is a vdW radius representing the distance at which the inter-atom potential is εij, and rij is the distance between a pair of atoms. For unlike atom pairs i and j, a 6th-order combining rule [26-29, 33] is used to calculate the off-diagonal parameters:

εij=2εiεjσi3σj3σi6+σj6,andσij=(σi6+σj62)16.
(10)

In the method used in previous publications regarding coarse graining [14-22], parameterization of the vdW interactions is conducted for every pair of CG beads in the system and pair-wise tabulated values of the potential are used in the calculation of non-bonded interactions. This results in the number of paired comparisons (Np) being

Np=i=1NCG(i)
(11)

where NCG is the number of CG beads in the system. This approach is satisfactory for polymers with relatively simple monomeric structures that can be represented with only a few different CG super-atoms. However, for more complex polymers, handling vdW parameterization in this manner rapidly becomes much more computationally expensive. For example, our CG model of poly(DTB succinate) contains five types of beads, and thus the number of combination of all pairs of beeds is fifteen. For more complex polymers in the TDPA family, there will be many more combinations and thus parameterization using the previously developed methods will be prohibitively expensive. Therefore, in this work, we developed an alternative, more efficient approach to parameterize the nonbonded interactions. We first generated vdW parameters for each type of bead in the CG model of poly(DTB succinate) by the following approach. Compounds that have similar chemical structure as the corresponding super-atom entities defined in the CG model (see Fig. 1) were first identified. The five types of low molecular weight compounds and their properties in the liquid state are listed in Table 2. Mono-phase systems for each of these five liquid-phase compounds were then built in a 25×25×25 Å3 periodic cube. For each liquid system, an all-atom NVT MD simulation was conducted at the corresponding liquid-state temperature given in Table 2 to generate a 2.0 ns trajectory.

The results from each MD simulation were then analyzed to generate a radial distribution function (RDF, or g(r)AA with r being the distance away from the center of mass of a reference molecule) for each model compound based on the center of mass of the molecules. These RDFs describe how, on average, the particles in a system radially pack around each other, thus providing an effective way to characterize their intermolecular interactions. These mono-phase systems were then coarse grained by replacing each molecule with its representative single CG bead, with the molecular weight and the center of mass of the CG bead and its atomistic counterpart being the same. An MD simulation was then performed for each of the five CG systems to generate another 2.0 ns trajectory, with the potential energy function used to determine the interactions between the CG beads in the system initially being defined as the potential of mean force, F(r), corresponding to the Boltzmann inversion of the g(r)AA determined from the all-atom MD simulation, or

F(r)=kBTlng(r)AA.
(12)

Simulating the interactions between the beads of a CG system with F(r) yields a corresponding RDF for the CG models, or g(r)CG, which can then be compared with the target RDF, or g(r)AA, to determine how well the vdW parameters CG model match the results of the all-atom simulations. Although F(r) is a free energy and not a potential energy (except for the case of zero density), it is usually sufficient to serve as the initial guess for the potential energy function of the CG models, or E0(r), from which, the parameters of σ0 and ε0 are estimated. Following the methods described by Soper [34], the CG parameters of σ and ε were optimized by an iterative Boltzmann inversion scheme, which uses the differences in the potentials of the mean force calculated by Eqn. (12) to iteratively improve the effective CG potential. The iterative scheme is accomplished by

Ei+1(r)=Ei(r)+kBTlngi(r)CGg(r)AA,
(13)

where Ei+1(r) and Ei(r) are CG potentials of successive iteration steps, which are used to replace the initial potential that was represented by F(r) from Eqn. (12). Valid solutions of σ and ε were obtained when the iterative procedure converged. A merit function was then used to measure the agreement between gi(r)CG and g(r)AA, which was defined by

f=ω(r)[g(r)AAgi(r)CG]2dr,
(14)

where ω (r) is a weighting function, defined as ω(r) = exp(-r).

Using this new approach, vdW parameters were first calculated for each of the five CG beads when they were interacting with similar bead types, and parameters for unlike pairs of beads were then calculated using the combining rule given by Eqn. (10).

2.2 Setup of Simulations

An atomistic simulation based on a hybrid scheme of MD and Monte Carlo (MC) was first performed for a DTB succinate tetramer to get the target distributions of bond length, bond angle and dihedral angle that were subsequently used to derive the CG parameters. Starting from an extended structure of DTB succinate tetramer, the hybrid MD/MC calculation was conducted in the NVT ensemble at 300 K with a 14 Å cut-off value being used for non-bonded interactions. Each MD/MC cycle was comprised of 1.0 ps MD followed by 1000 steps of MC. The MD simulation was conducted with the velocity Verlet algorithm [35] and a 0.25 fs time step, while the dihedral angles were randomly sampled by the MC procedure. In all, 2 × 104 cycles of MD/MC were conducted, corresponding to 5.0 ns MD and 2 × 107 samples of MC, respectively. The final structure of the atomistic MC/MD simulation was then saved and mapped onto a CG structure. Each of these structures (i.e., the all-atom and the CG model) was then used in a regular MD simulation to generate a 12 ns trajectory for each representation of the tetramer. The resulting trajectories were then used to validate and assess the degree of performance enhancement provided by the CG model by calculating and comparing the radius of gyration and the diffusivity of the CG representation with the corresponding properties of the atomistic representation of the tetramer.

Once the parameterization of the CG model was validated, a CG model of DTB succinate was further used in an MD simulation to relax a dense-phase system composed of eight chains of poly(DTB succinate), with each polymer chain containing 50 repeat units. To prepare the starting configuration, an MD calculation was performed at 600 K to relax the chains, which were initially structured in their fully extended conformations. A 10 kcal/mol harmonic potential was imposed to compress and then constrain all chain segments within a cubic box with dimensions 64.43 × 64.43 × 64.43 Å3. These dimensions were set to provide a system density of about 1.2 g/cm3, which is similar to that of bulk poly(DTB succinate) at 300 K temperature [36]. This constrained high temperature CG polymer model was then used as the initial structure for an MD simulation conducted at 300K. The polymers were first energy-minimized for 1000 steps using the adopted-basis Newton-Raphson (ABNR) method [25], heated to 300 K in 25 ps, and then thermally equilibrated for 50 ps. A 10 ns trajectory was then generated in the production phase of the MD simulation using 2 fs integral time step to generate the final relaxed, dense-phase CG model of the amorphous polymer structure.

2.3 Reverse Mapping from the CG Model to the All-Atom Dense-Phase Polymer Model

The resulting CG structure from the 10 ns MD simulation was used to rebuild the atomistic structure using a reverse-mapping procedure. In the present implementation of this procedure, the atomistic structure was rebuilt for each chain by sequences of successive superposition between adjacent CG beads starting from one end-group of each polymer chain. The atomistic counterparts for each bead were superimposed onto the corresponding CG bead by overlapping their centers of mass, with the orientation of the all-atom model generally set to correspond to the starting and ending atom for each chain segment. To form a chemical bond between a superimposed and a newly inserted atomistic unit, the inserted superimposed atomistic unit was allowed to rotate freely but with fixed center of mass and with a 10 kcal/mol attractive harmonic constraint being imposed between the two atoms that needed to be connected to join the adjacent CG elements. ABNR energy-minimization was then conducted until the distance between the two atoms became less than 1.5 ´. A new bond was then formed between the adjacent atom pair, following which a further 100 steps of energy minimization were performed to optimize the new structure. At the end of this reverse-mapping procedure, 500 additional steps of energy minimization were then conducted to optimize the global structure of the polymer, thus finalizing the development of the all-atom, dense phase model of the polymer.

3. RESULTS AND DISCUSSION

3.1 Coarse Graining

The coarse graining procedure that was used for the development of the CG parameters to represent the polymer chain was performed based on a single chain of DTB succinate tetramer. The hybrid MD/MC scheme that was applied to initially sample the atomistic behavior of this oligmer resulted in an average MC sampling acceptance ratio of about 40% and provided an efficient approach to equilibrate the molecule. The resulting evolution of the potential energy of this system as a function of MC/MD sampling cycles is plotted in Fig. 2. As shown, the energy of the system is equilibrated after only a few MC/MD cycles and then fluctuates about 7.0 ± 17.7 kcal/mol throughout the simulation.

Fig. 2
Potential energy (E) vs. number of MD/MC sampling cycles (Ncycle) in the atomistic simulation of a DTB succinate tetramer.

Based on the last 104 cycles of the atomistic simulation, the distributions of structural factors were determined, from which the bonded parameters of the CG force field were derived. The six distributions of distances between adjacent super-atoms obtained from atomistic calculation (scattered points) and double Gaussian function fittings (lines) are shown in Fig. 3. The fitted parameters for bond stretching interaction [Eqn. (5)] are given in Table 3. The eight distributions of bond angles between three successive super-atoms obtained from the atomistic calculation (scattered points) and double Gaussian function fittings (lines) are shown in Fig. 4. The fitted parameters for the bending interactions [Eqn. (7)] are given in Table 4. As clearly evident from these figures, the double Gaussian function is effective in fitting all of the distributions of bond length and angle.

Fig. 3
Histogram of the length between adjacent CG beads obtained from the atomistic simulation (scattered points) and the nonlinear fitting based on double Gaussian functions (lines). The bond types are defined in Table 1.
Fig. 4
Histogram of the angle between three successive CG beads obtained from the atomistic simulation (scattered points) and the nonlinear fitting based on double Gaussian functions (lines). The bond angles are defined in Table 1.
Table 3
The bond stretching interaction parameters for the CG model of DTB succinate fitted by the double Gaussian function [Eqn. (4)] for the atomistic bond distributions.
Table 4
The bending interaction parameters for the CG model of DTB succinate fitted by the double Gaussian function for the atomistic bond angle distributions.

For the five types of beads (Fig. 1) in the CG model of poly(DTB succinate), the vdW parameters in Eqn. (9) were obtained by mapping the radial distribution function, RDF, of all-atom models of mono-dispersed compounds (Table 2) in their liquid state that have similar chemical structure as the super-atoms in the polymer. Mathematically, to calculate RDFs, a molecule in the system is chosen as a reference and a series of concentric spheres are drawn around it. Considering a spherical shell of thickness δr at a distance r from the center of mass of the reference molecule, the probability of finding another molecule is calculated by

g(r)=1dn(r)4πr2δr,
(15)

where g(r) represents the RDF, n(r) is the mean number of molecules in the shell (counted based on the position of the center of mass), and d is the mean bulk density. The results of the calculated values of g(r) for each of the compounds shown in Table 2 are plotted in Fig. 5 with filled circles being connected by dotted lines. Within a zone of short separations (i.e., small r), g(r) is zero, indicating the effective radii of the molecules where they sterically resist overlap with one another. The occurrence of peaks at long range indicates that the molecules pack around each other in shells of neighbors. At long range, g(r) tends to a value of 1.0, thus describing the average density at this range. Starting from the initial RDF obtained from each liquid compound, the nonbonded interaction parameters were calculated by the iterative Boltzmann inversion scheme defined by Eqns. (12) and (13). Consistently with published work [20, 31], all RDFs converged after about four iterations. The RDFs at the initial and final stages are plotted in Fig. 5 with dashed and solid lines, respectively. The best fit result for each type of CG bead associated with the value of the merit function [Eq. (14)] when the RDF converged are summarized in Table 5. Although reproducing the RDFs of the all-atom models with the CG models was clearly more difficult than reproducing the behavior of the bonded interactions (as expected), the RDFs from the CG simulations were able to be adjusted to closely match the position and height of the main peaks at small separation for each of the five representative model systems. The RDFs of MeF and PhF beads were particularly hard to be represented by their respective CG models. These deviations are caused by the asymmetric geometry of the molecule. For example, the atomistic structure of PhF is basically planar due to the presence of the phenyl ring. The strong hydrophobic interaction between nonbonded phenyl rings leads to close packing in parallel orientation, resulting in a sharp peak in the RDF at about 2.2 ´ of separation. Broad peaks then also occur at about 5.5 ´ distance, corresponding to the length of the molecule. In the current mapping scheme, in which all super-atoms were described as spherical beads, the vdW radius of the PhF bead was approximated by reproducing the first peak of the RDF. As presented below, comparisons of the radius of gyration between the all-atom and the CG models of the DTB succinate tetramer show that the developed CG model is able to do an excellent job in representing the overall structural behavior of the polymer chain, thus indicating that the values of CG parameters that were generated by these methods were reasonable.

Fig. 5
Optimization of the intermolecular RDF for the five liquid compounds shown in Table 2. The initial guess of RDF is shown by the dashed line. The quality of the trial RDFs improves following successive iterations, and the final RDF (solid line) matche ...
Table 5
The optimized CG parameters for non-bonded interaction described by Lennard-Jones 6-9 potential [Eqn. (9)]. For each type of bead, the value of the merit function, f [Eqn. (14)] is presented for when the radial distribution function (RDF) was converged. ...

Lastly, the target distributions of the nine types of dihedral angles in DTB succinate were calculated from the atomistic model of DTB succinate tetramer. The results from these calculations are shown in Fig. 6 with filled circles. The force constants, Ki, and phase angles, δi (i = 1, 3), were optimized so that, based on the torsional potential defined by Eqn. (8), the distributions of dihedral angle calculated by the CG model could reproduce the target distributions. The converged CG distributions are plotted in Fig. 6 with solid lines and the optimized parameters corresponding to each angle are given in Table 6.

Fig. 6
Distribution of the CG dihedral angles obtained from the atomistic simulation (scattered points) and the nonlinear fitting of the potential based on Eqn. (8) (lines). The dihedral angle types are defined in Table 1.
Table 6
The optimized CG parameters for dihedral interaction term described by Eqn. (8).

3.2 Mapping Results of Static and Dynamic Properties

The development of a reasonable CG model requires that it not only reproduces the local properties of the polymer, but it also reproduces the global static properties of the polymer. Of course, since the primary purpose of developing a CG model is to improve computational efficiency, it is also important to assess just how much efficiency is gained by implementing the CG scheme. To address these issues, the root-mean-square radius of gyration was calculated for atomistic DTB succinate tetramer from the MD/MC calculation. Furthermore, separate MD simulations were conducted using the atomistic and the CG models of DTB succinate tetramer. Each simulation generated a 12 ns trajectory, based on which the root-mean square radius of gyration of the CG model and the mean-square displacement of the center of mass of the tetramer represented by both atomistic and CG models were extracted.

The root—mean—square radius of gyration (<RG2>1/2) is a measure of the size of a polymer chain, and is defined as the root—mean—square distance of atoms from their common center of gravity. Mathematically, <RG2>1/2 is calculated by,

RG212=1Ni=1N(RiRG)212,andRG=1Ni=1NRi,
(16)

where N is the chain length, Ri and RG are the positions of an atom and the center of gravity of the polymer, respectively. The calculated values of <RG2> from the atomistic and CG simulations are given in Table 7. Note that the CG data were calculated with two sets of parameters of PhF bead, i.e. (1) σ = 2.22 ´ and ε = 0.4775 kcal/mol, and (2) σ = 5.43 ´ and ε = 0.1095 kcal/mol. Parameter set (1) emphasizes the short-range peak in the RDF while set (2) emphasizes the broad long-range peak in the RDF of PhF (Fig. 5). Table 7 shows that the CG result with PhF parameter set (1) agrees more closely sizes with the atomistic result, whereas PhF parameter set (2) leads to larger deviation from the atomistic data. Therefore, PhF parameter set (1) was adopted in the development of the CG model for poly(DTB succinate).

Table 7
The radius of gyration of DTB succinate tetramer calculated using atomistic and CG models, respectively.

A useful measure to evaluate the improvement of the efficiency for the CG model is to compare the molecular diffusivity obtained from the atomistic and CG models, respectively. According to the treatment in polymer dynamics [37], a flexible chain at equilibrium state can be approximated as a random coil. The trajectory of the chain is equivalent to that of a Brownian particle; hence the diffusivity D can be defined in terms of mean-square displacement of the center of mass of the coil by

D=limt16t[R(t)RG(0)]2,
(17)

where t is time, RG(0) and RG(t) are the centers of mass at time zero and t, respectively. We estimated D from MD simulations for both the atomistic and the CG models by following the center of mass of the tetramer versus time, with the results shown in Fig. 7. D is estimated directly from the slope of each profile, and the results are given in Table 8. As shown, the diffusivity of the CG model scales about 120 times larger than that of the atomistic model, thus indicating that the CG model provides a over two orders of magnitude in speedup compared to the corresponding atomistic model.

Fig. 7
Displacement of the center of mass of DTB succinate tetramer obtained from atomistic and CG models, respectively.
Table 8
The results of diffusivity (D) of DTB succinate tetramer calculated based on atomistic and CG models, respectively.

3.3 Reverse Mapping

As the final component of this program, the developed CG model was applied to generate a relaxed, dense-phase all-atom amorphous polymer model of poly(DTB succinate) as represented by a set of eight polymer chains, each with 50 monomeric repeat units. The CG polymer model was first relaxed by an MD simulation at 600 K, energy minimized, and then re-equilibrated at 300 K. A final 10 ns MD was then conducted at 300 K to further relax the polymer chains. The atomistic structure of the system was then reconstructed by reverse mapping from the final CG structure after the 10 ns MD simulation for system relaxation. Fig. 8A shows one of the CG chains and its atomistic counterpart. The whole system contained in a cubic box is shown in Fig. 8B, with an overall density of 1.2 g/cm3, which closely matches the experimentally measured density of poly(DTB succinate) [36].

Fig. 8Fig. 8
(A) A CG model of poly(DTB succinate) and the corresponding atomistic model obtained from the reverse mapping scheme. (B) The bulk structure of eight chains (shown in different colors) of 50-mer poly(DTB succinate) rebuilt from a CG structure. The density ...

These results demonstrate that the developed CG methods can be used to construct an all-atom dense-phase model of a complex branched amorphous polymer starting from a given all-atom description of the monomer unit of the polymer chain. It must be recognized, however, that it is unlikely that the resulting structure, as shown in Fig. 8B, represents a completely equilibrated model of the amorphous polymer, with the implementation of additional advanced sampling algorithms being needed for final system equilibration. While this is acknowledged, the dense-phase models that can be generated by the presented CG methods provide an important starting point from which subsequent advanced sampling methods can be applied, and we are currently working the development and application of advanced sampling methods for this purpose.

4. CONCLUSION

While all-atom molecular models of single oligomeric chains of a given amorphous polymer can be relatively easily constructed, it is a much greater challenge to develop polymer models that realistically represent the dense-phase structure of an amorphous polymer such that they can then be used to simulate polymer behavior for various areas of application. To address this challenge, we have developed an approach to generate a CG model of a polymer chain that enables the general structural characteristics of the polymer chain to be represented, but in a manner that greatly reduces the number of degrees of freedom in the system. This enables the CG polymer model to be readily used to generate a relaxed, dense-phase CG representation of the bulk polymer in an efficient manner, with reverse-mapping and energy minimization then applied to translate the CG model back into an all-atom representation of the system to finally provide a representative relaxed, dense-phase, all-atom model of the amorphous polymer.

In this paper, we demonstrate the application of these methods for a complex branched amorphous polymer, poly(DTB succinate). Beginning with an all-atom model of a DTB succinate oligmer, these procedures enabled us to construct a relaxed, dense-phase, all-atom molecular model of this polymer with a density matching the experimentally measured density of the bulk polymer at 300 K temperature. It is recognized that these initial dense-phase models do not yet represent a truly equilibrated structure of the polymer. Additional advanced sampling algorithms are needed and are being planned to attain this level of equilibration, with the equilibrated polymer structures to be validated by comparing RDFs from the polymer model with RDFs obtained by wide-angle x-ray scattering (WAXS) measurements [38, 39]. The present work sets the stage from which these further studies can be implemented to attain our final objective of producing accurate dense-phase all-atom models of complex amorphous polymers, which can then be used as a guide for polymer design to optimize performance for targeted applications.

ACKNOWLEDGEMENT

Financial support for this work was provided by RESBIO, a Research Resource funded by the National Institutes of health under NIH grant EB001046, the New Jersey Center for Biomaterials, and Rutgers University. The authors also acknowledge the RESBIO program at Rutgers University for providing the experimental density of poly(DTB succinate).

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

REFERENCES

(1) Xue L, Greisler HP. J. Vasc. Surg. 2003;37:472. [PubMed]
(2) Panyam J, Labhasetwar V. Adv. Drug Delivery Rev. 2003;55:329. [PubMed]
(3) Langer R, Vacanti J. Science. 1995;260:920. [PubMed]
(4) Kumar N, Ravikumar MN, Domb AJ. Adv. Drug Delivery Rev. 2001;53:23. [PubMed]
(5) Brocchini J, James K, Tangpasuthadol V, Kohn J. J Biomed Mater Res. 1998;42:66. [PubMed]
(6) Brocchini S, James K, Tangpasuthadol V. J Am Chem Soc. 1997;119:4553.
(7) Latour RA. Biomaterials: Protein-Surface Interactions. In: Wnek GE, Bowlin GL, editors. The Encyclopedia of Biomaterials and Bioengineering. 2nd Edition Vol. 1. Informa Healthcare; 2008. p. 270.
(8) Kholodovych V, Smith JR, Knight D, Abramson S, Kohn J, Welsh WJ. Polymer. 2004;45:7367.
(9) Smith JR, Seyda A, Weber N, Knight D, Abramson S, Kohn J. Macromol Rapid Commun. 2004;25:127.
(10) Sousa A, Aitouchen A, Libera M. Ultramicroscopy. 2006;106:130. [PubMed]
(11) Sharma RI, Kohn J, Moghe PV. J Biomed Mater Res. 2004;69:114. [PubMed]
(12) Gubskaya AV, Kholodovych V, Knight D, Kohn J, Welsh WJ. Polymer. 2007;48:5788. [PMC free article] [PubMed]
(13) Baumgartner A. Annu Rev Phys Chem. 1984;35:419.
(14) Shelley JC, Shelley MY, Reeder RC, Bandyopadhyay S, Klein ML. J Phys Chem B. 2001;105:4464.
(15) Shelley JC, Shelley MY, Reeder RC, Bandyopadhyay S, Moore PB, Klein ML. J Phys Chem B. 2001;105:9785.
(16) Guerrault X, Rousseau B, Farago J. J Chem Phys. 2004;121:6538. [PubMed]
(17) Tschöp W, Kremer K, Batoulis J, Bürger T, Hahn O. Acta Polymer. 1998;49:61.
(18) León S, van der Vegt N, Delle Site L, Kremer K. Macromolecules. 2005;38:8078.
(19) Hess B, León S, van der Vegt N, Kremer K. Soft Matter. 2006;2:409.
(20) Li X, Kou D, Rao S, Liang H. J Chem Phys. 2006;124:204909. [PubMed]
(21) Srinivas G, Discher DE, Klein ML. Nature Materials. 2003;3:638. [PubMed]
(22) Ortiz V, Nielsen SO, Klein ML, Discher DE. J Polym Sci B: Polym Phys. 2006;44:1907.
(23) Murat M, Kremer K. J Chem Phys. 1998;108:4340.
(24) Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. J Comput Chem. 1983;4:187.
(25) MacKerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher WE, Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M. J Phys Chem B. 1998;102:3586. [PubMed]
(26) Hariharan PC, Pople JA. Theor Chim Acta. 1973;28:213.
(27) Francl MM, Pietro WJ, Hehre WJ, Binkley JS, Gordon MS, DeFrees DJ, Pople JA. J Chem Phys. 1982;77:3654.
(28) Dinur U, Hagler AT. In: Reviews in Computational Chemistry. Lipkowitz KB, Boyd DB, editors. Vol. 2. VCH Publishers; New York: 1991. p. 99.
(29) Maple JR, Hwang MJ, Stockfisch TP, Dinur U, Waldman M, Ewig CS, Hagler AT. J Comp Chem. 1994;15:162.
(30) Müller-Plathe F. Chem Phys Chem. 2002;3:754.
(31) Reith D, Meyer H, Müller-Plathe F. Macromolecules. 2001;34:2335.
(32) Tschöp W, Kremer K, Hahn O, Batoulis J, Bürger J. Acta Polymer. 1998;49:61.
(33) Waldman M, Hagler AT. J Comput Chem. 1993;14:1077.
(34) Soper AK. Chem Phys. 1996;202:295.
(35) Allen MP, Tildesley DJ. Computer Simiulation of Liquids. Clarendon Press; Oxford: p. 1987.
(36) Internal correspondence from the group of Dr. Knight D, Rutgers University.
(37) Doi M, Edwards SF. The theory of polymer dynamics. Clarendon; Oxford: 1986.
(38) Habenschuss A, Tsige M, Curro JG, Grest GS, Nath SK. Macromolecules. 2007;40:7036.
(39) Mitchell GR, Rosi-Schwartz B, Ward DJ, Warner M. Local Order in Polymer Glasses and Melts. Philosophical Transactions: Physical Sciences and Engineering. 1994;348:97.