Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2742384

Formats

Article sections

Authors

Related links

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.055PMCID: PMC2742384

NIHMSID: NIHMS128961

Department of Bioengineering, Clemson University, Clemson, South Carolina 29634

See other articles in PMC that cite the published article.

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 10^{2}, 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.

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 *N*^{3}, leading to more than *N*^{4} 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.

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.

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 (Et_{1} and Et_{2}), 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 Et_{1} and Et_{2} 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.

(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 **...**

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={E}_{\text{bonded}}+{E}_{\text{non-bonded}},$$

(1)

where *E*_{bonded} and *E*_{non-bonded} are bonded and nonbonded potential terms, respectively, which are decomposed as,

$$\begin{array}{c}{E}_{\text{bonded}}={E}_{\text{stretch}}+{E}_{\text{bend}}+{E}_{\text{torsion}}+{E}_{\text{cross}}\hfill \\ {E}_{\text{non-bonded}}={E}_{\mathit{vdW}}+{E}_{\text{Coulomb}}\hfill \end{array},$$

(2)

where *E*_{stretch}, *E*_{bend} and *E*_{torsion} are the bond-length-stretching, angle-bending, and bond-rotation energies, respectively; the cross term, *E*_{cross}, accounts for the coupling between individual bonded interaction; *E*_{vdW} accounts for the excluded volume repulsive as well as the intermolecular attractive forces between non-bonded atoms; and *E*_{Coulomb} is the electrostatic potential. The CG force field is developed in a similar manner with a few modifications. Specifically, the *E*_{cross} 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 *E*_{stretch}, *E*_{bend}, *E*_{torsion}, and *E*_{vdW}. 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 *E*_{stretch} → *E*_{bend} → *E*_{vdW} → *E*_{torsion}. 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.

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

$$P\left(l\right)={Z}_{l}p\left(l\right)\u2215{l}^{2}=\mathrm{exp}[-{E}_{\text{stretch}}\u2215{k}_{B}T],$$

(3)

where *l* is the distance between successive beads, Z_{l} is a normalizing factor, *l*^{2} 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\left(l\right)=\frac{{A}_{1}}{{w}_{1}}\mathrm{exp}\left(\frac{-2{(l-{l}_{c1})}^{2}}{{w}_{1}^{2}}\right)+\frac{{A}_{2}}{{w}_{2}}\mathrm{exp}\left(\frac{-2{(l-{l}_{c2})}^{2}}{{w}_{2}^{2}}\right)=\mathrm{exp}[-{E}_{\text{stretch}}\u2215{k}_{B}T],$$

(4)

where *A*_{1}, *A*_{2}, *l _{c}*

$${E}_{\text{stretch}}={k}_{B}T\left\{{K}_{b1}{(l-{l}_{c1})}^{2}-\mathrm{ln}[1+{K}_{b3}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}[{K}_{b1}{(l-{l}_{c1})}^{2}-{K}_{b2}{(l-{l}_{c2})}^{2}\left]\right]\right\},$$

(5)

where the energetic parameters are defined as: *K _{b}*

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\left(\theta \right)={Z}_{n}p\left(\theta \right)\u2215\mathrm{sin}\left(\theta \right)=\mathrm{exp}\left[-{E}_{\text{bend}}\u2215{k}_{B}T\right],$$

(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 *E _{bend}*; i.e.,

$${E}_{\text{bend}}={k}_{B}T\left\{{K}_{\theta 1}{(\theta -{\theta}_{c1})}^{2}-\mathrm{ln}[1+{K}_{\theta 3}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}[{K}_{1}{(\theta -{\theta}_{c1})}^{2}-{K}_{\theta 2}{(\theta -{\theta}_{c2})}^{2}\left]\right]\right\},$$

(7)

where *K _{θ}*

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

$${E}_{\text{torsion}}\left(\varphi \right)=\sum _{n=1}^{3}{K}_{n}\left[1+\mathrm{cos}(n\varphi -{\delta}_{n})\right],$$

(8)

where *ϕ* is the dihedral angle, *K _{i}* and

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.,

$${E}_{\mathit{vdw}}={\epsilon}_{ij}\left[2{\left(\frac{{\sigma}_{ij}}{{r}_{ij}}\right)}^{9}-3{\left(\frac{{\sigma}_{ij}}{{r}_{ij}}\right)}^{6}\right],$$

(9)

where *ε _{ij}* is the depth of the potential well,

$${\epsilon}_{ij}=2\sqrt{{\epsilon}_{i}{\epsilon}_{j}}\frac{{\sigma}_{i}^{3}{\sigma}_{j}^{3}}{{\sigma}_{i}^{6}+{\sigma}_{j}^{6}},\phantom{\rule{thickmathspace}{0ex}}\text{and}\phantom{\rule{thickmathspace}{0ex}}{\sigma}_{ij}={\left(\frac{{\sigma}_{i}^{6}+{\sigma}_{j}^{6}}{2}\right)}^{1\u22156}.$$

(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 (*N _{p}*) being

$${N}_{p}=\sum _{i=1}^{{N}_{\mathit{CG}}}\left(i\right)$$

(11)

where *N _{CG}* 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 Å

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\left(r\right)=-{k}_{B}T\mathrm{ln}\phantom{\rule{thinmathspace}{0ex}}g(r{)}_{\mathit{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 *E*_{0}(*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

$${E}_{i+1}\left(r\right)={E}_{i}\left(r\right)+{k}_{B}T\mathrm{ln}\frac{{g}_{i}(r{)}_{\mathit{CG}}}{g(r{)}_{\mathit{AA}}},$$

(13)

where *E _{i}*

$$f=\int \omega \left(r\right){\left[g(r{)}_{\mathit{AA}}-{g}_{i}(r{)}_{\mathit{CG}}\right]}^{2}dr,$$

(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).

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 × 10^{4} cycles of MD/MC were conducted, corresponding to 5.0 ns MD and 2 × 10^{7} 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/cm^{3}, 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.

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 $\stackrel{\xb4}{\mathrm{\u212b}}$. 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.

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.

Potential energy (*E*) vs. number of MD/MC sampling cycles (N_{cycle}) in the atomistic simulation of a DTB succinate tetramer.

Based on the last 10^{4} 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.

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\left(r\right)=\frac{1}{d}\frac{n\left(r\right)}{4\pi {r}^{2}\delta 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 $\stackrel{\xb4}{\mathrm{\u212b}}$ of separation. Broad peaks then also occur at about 5.5 $\stackrel{\xb4}{\mathrm{\u212b}}$ 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.

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, *K _{i}*, and phase angles,

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 (<*R _{G}*

$${\langle {R}_{G}^{2}\rangle}^{1\u22152}=\frac{1}{N}\sum _{i=1}^{N}{\langle {({\overrightarrow{R}}_{i}-{\overrightarrow{R}}_{G})}^{2}\rangle}^{1\u22152},\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{\overrightarrow{R}}_{G}=\frac{1}{N}\sum _{i=1}^{N}{\overrightarrow{R}}_{i},$$

(16)

where *N* is the chain length, ${\overrightarrow{R}}_{i}$ and ${\overrightarrow{R}}_{G}$ are the positions of an atom and the center of gravity of the polymer, respectively. The calculated values of <*R*_{G}^{2}> 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 $\stackrel{\xb4}{\mathrm{\u212b}}$ and *ε* = 0.4775 kcal/mol, and (2) *σ* = 5.43 $\stackrel{\xb4}{\mathrm{\u212b}}$ 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).

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=\underset{\mathrm{t}\to \infty}{\mathrm{lim}}\frac{1}{6t}\langle {\left[\overrightarrow{R}\left(t\right)-{\overrightarrow{R}}_{G}\left(0\right)\right]}^{2}\rangle ,$$

(17)

where *t* is time, ${\overrightarrow{R}}_{G}\left(0\right)$ and ${\overrightarrow{R}}_{G}\left(t\right)$ 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.

Displacement of the center of mass of DTB succinate tetramer obtained from atomistic and CG models, respectively.

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/cm^{3}, which closely matches the experimentally measured density of poly(DTB succinate) [36].

(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.

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.

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).

**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.

(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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |