Exploring the Role of Hydroxy- and Phosphate-Terminated <i>cis</i>-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber: Insights from Molecular Dynamics Simulations.

Dixit M; Taniguchi T

doi:10.1021/acspolymersau.4c00019

Exploring the Role of Hydroxy- and Phosphate-Terminated cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber: Insights from Molecular Dynamics Simulations.

Dixit M ¹,

Taniguchi T ¹

Affiliations

1. Graduate School of Engineering, Kyoto University, Nishikyo-ku, Kyoto 615-8510, Japan.
Authors
Dixit M¹
Taniguchi T¹
(2 authors)

ORCIDs linked to this article

ACS Polymers Au, 15 May 2024, 4(4):273-288
https://doi.org/10.1021/acspolymersau.4c00019 PMID: 39156555 PMCID: PMC11328332

This article is in the Europe PMC Open access subset. Refer to the copyright information in the article for licensing details.

Free full text in Europe PMC

Abstract

This study elucidates the pivotal role of terminal structures in cis-1,4-polyisoprene (PI) chains, contributing to the exceptional mechanical properties of Hevea natural rubber (NR). NR's unique networking structure, crucial for crack resistance, elasticity, and strain-induced crystallization, involves two terminal groups, ω and α. The proposed ω terminal structure is dimethyl allyl-(trans-1,4-isoprene)₂, and α terminals exist in various forms, including hydroxy, ester, and phosphate groups. Among others, we investigated three types of cis-1,4-PI with different terminal combinations: _HPI_H (pure PI with H terminal), _ωPI_α6 (PI with ω and α6 terminals), and _ωPI_PO₄ (PI with ω and PO₄ terminals) and revealed significant dynamics variations. Hydrogen bonds between α6 and α6 and PO₄ and PO₄ residues in _ωPI_α6 and _ωPI_PO₄ systems induce slower dynamics of hydroxy- and phosphate-terminated PI chains. Associations between α6 and α6 and PO₄ and PO₄ terminals are markedly stronger than ω and ω, and hydrogen terminals in _HPI_H and _ω PI_α6,PO₄ systems. Phosphate terminals exhibit a stronger mutual association than hydroxy terminals. Potentials of mean force analysis and cluster-formation-fraction computations reveal stable clusters in _ωPI_α6 and _ωPI_PO₄ , supporting the formation of polar aggregates (physical junction points). Notably, phosphate terminal groups facilitate large and highly stable phosphate polar aggregates, crucial for the natural networking structure responsible for NR's outstanding mechanical properties compared to synthetic PI rubber. This comprehensive investigation provides valuable insights into the role of terminal groups in cis-1,4-PI melt systems and their profound impact on the mechanical properties of NR.

Free full text

ACS Polym Au. 2024 Aug 14; 4(4): 273–288.

Published online 2024 May 15. https://doi.org/10.1021/acspolymersau.4c00019

PMCID: PMC11328332

PMID: 39156555

Exploring the Role of Hydroxy- and Phosphate-Terminated cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber: Insights from Molecular Dynamics Simulations

Mayank Dixit* and Takashi Taniguchi*

Author information Article notes Copyright and License information Disclaimer

Associated Data

Supplementary Materials: lg4c00019_si_001.pdf
lg4c00019_si_001.pdf (14M)
lg4c00019_si_002.mp4
lg4c00019_si_002.mp4 (52M)
lg4c00019_si_003.mp4
lg4c00019_si_003.mp4 (51M)

Abstract

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0013.jpg

This study elucidates the pivotal role of terminal structures in cis-1,4-polyisoprene (PI) chains, contributing to the exceptional mechanical properties of Hevea natural rubber (NR). NR’s unique networking structure, crucial for crack resistance, elasticity, and strain-induced crystallization, involves two terminal groups, ω and α. The proposed ω terminal structure is dimethyl allyl-(trans-1,4-isoprene)₂, and α terminals exist in various forms, including hydroxy, ester, and phosphate groups. Among others, we investigated three types of cis-1,4-PI with different terminal combinations: _HPI_H (pure PI with H terminal), _ωPI_α6 (PI with ω and α6 terminals), and _ωPI_PO₄ (PI with ω and PO₄ terminals) and revealed significant dynamics variations. Hydrogen bonds between α6 and α6 and PO₄ and PO₄ residues in _ωPI_α6 and _ωPI_PO₄ systems induce slower dynamics of hydroxy- and phosphate-terminated PI chains. Associations between α6 and α6 and PO₄ and PO₄ terminals are markedly stronger than ω and ω, and hydrogen terminals in _HPI_H and _ωPI_α6,PO₄ systems. Phosphate terminals exhibit a stronger mutual association than hydroxy terminals. Potentials of mean force analysis and cluster-formation-fraction computations reveal stable clusters in _ωPI_α6 and _ωPI_PO₄, supporting the formation of polar aggregates (physical junction points). Notably, phosphate terminal groups facilitate large and highly stable phosphate polar aggregates, crucial for the natural networking structure responsible for NR’s outstanding mechanical properties compared to synthetic PI rubber. This comprehensive investigation provides valuable insights into the role of terminal groups in cis-1,4-PI melt systems and their profound impact on the mechanical properties of NR.

Keywords: phosphate terminal group, large size cluster, potentials of mean force, hydrogen bond, polar aggregates

1. Introduction

The unique networking structure in natural rubber (NR) imparts superior chemical and mechanical properties compared to synthetic cis-1,4-polyisoprene rubber (PI).¹⁻¹² Although synthetic PI rubber has been adopted as an alternative due to limited NR production, its comprehensive properties still lag behind NR.¹³ The exceptional properties of NR over PI rubber arise from the natural networking structure formed between nonrubber components and the terminal groups of cis-1,4-PI chains.^3,11,12 Therefore, elucidating the role of terminal groups in NR can provide valuable insights for designing new rubber materials. Hevea NR consists of rubber particles, protein, lipids, carbohydrates, metal ions, and other components.^2,14−19 The main constituents of Hevea NR are cis-1,4-PI and two types of terminal groups.^{2,14,15,20,21} These terminal groups are referred to as ω-terminal and α-terminal. The ω-terminal consists of a dimethylallyl group and 2–3 trans isoprene. The α-terminals are considered to be categorized into those containing hydroxy, ester, and phosphate groups.^{1,3,16,20−25} Previous investigations have conducted structural analysis using solution NMR on the terminal units of NR subsequent to chemical treatment involving alterations such as deproteinization with enzymes and other chemicals.^2,26,27 These studies revealed the presence of monophosphate, diphosphate, and phospholipids at the α-terminus.^2,26−28 In contrast, Oouchi et al. devised a NMR-based approach to analyze terminal unit structures in commercial NR, either without chemical treatments or with mild treatments like acetone extraction to remove impurities. Their novel methods, integrating 2D-NMR, effectively suppressed signals from low-molecular-weight impurities. Employing these techniques, they detected NMR signals of terminal units in chemically untreated NR. Analysis of eight commercial H-NR samples consistently revealed six α-terminus terminating-end units, with none forming a phosphate ester.¹ Hence, there is a contradiction between these two experimental research studies. This disparity between experimental findings necessitates a thorough investigation into the molecular characteristics of these α-terminal groups in NR. It is also reported in the literature^2,3,24,29 that, in a melt of natural cis-1,4-PI, a physical network structure is formed through hydrogen bonds (HBs) in α-terminal and through hydrogen and ionic bonds with phospholipids, while ω-terminal bonds are formed with protein molecules. These branch points in the network play a vital role in polymer chain elongation during deformation. Previous research utilizing NMR has suggested that lipid molecules form physical junction points with α-terminal through hydrogen and ionic bonds.^{2,14,15,22,23,25} Recent experimental studies have also indicated the significance of physical junction points by terminal groups in the strain-induced crystallization (SIC) of NR.^18,30−32 Moreover, a synthesized phosphate-terminated PI rubber has shown improved mechanical properties compared to pure PI rubber due to multiple physical junction points mediated by phosphate terminal groups.³³ Recent findings demonstrate that phosphorylcholine groups function as cross-linking sites (physical junction points) in both unvulcanized and vulcanized rubber, thereby enhancing mechanical characteristics through the facilitation of SIC. Conversely, polymers featuring pendant hydroxyl groups exhibit sequence-specific SIC behaviors attributed to their nonaggregating nature.³⁴ However, the structural mechanism of physical junction point formation remains unclear. To address this knowledge gap, we perform coarse-grained and all-atom molecular dynamics (MD) simulations, investigating the role of ω-terminal and two types of α′-terminals (i.e., α6-terminal group and IP–PO₄ terminal group) in the formation of physical junction points by PI chains. In this work, we used the terminology α′-terminal instead of α-terminal because the terminal-containing PO₄ has not been confirmed yet by an experiment to be a terminal group of cis-1,4-PI taken from a rubber tree. Our multiscale approach is based on our recent work,²²⁻²⁵ where we developed a model for cis-1,4-PI chains with six types of terminal groups (i.e., from α1 to α6, according to their notation).¹ The α6-terminal group represents hydroxy-terminated 1,4-cis-isoprene [−CH₂–C(CH₃)=CH– CH₂–OH] and IP–PO₄] terminal group represents phosphate-terminated 1,4-cis-isoprene [−CH₂–C(CH₃)=CH–CH₂–H₂PO₄]. Therefore, there is structural similarity of monomer units in between α6 and IP–PO₄ terminal groups. Henceforth, the notation for the −H₂PO₄ group shall be simplified to −PO₄. Table 1 provides details of the three distinct melt systems with different terminal groups (see Figure Figure11). Molecular modeling and MD simulations have been invaluable in investigating the mechanical and chemical properties of polymers, complementing experimental studies.^24,35−62 In our recent research,²²⁻²⁵ we employed a multiscale approach to explore the interactions of terminal groups in the cis-1,4-PI chain, which was recently generated as reported by Choi et al.⁶³ In the current study, we focus on cis-1,4-PI with ω-terminal and two distinct types of α′ terminals, namely, the α6-terminal group and the PO₄ terminal group. In addition, as a reference melt system, we also consider a melt system of hydrogen-terminated cis-1,4-PI (_HPI_H). Namely, we investigate three different melt systems (_HPI_H), _ωPI_α6, and _ωPI_PO₄, each characterized by unique combinations of ω-terminal and α′-terminal groups (illustrated in Figure Figure11). Our approach builds upon the same multiscale model utilized in our previous work.²²⁻²⁵ Detailed information about each melt system is provided in Table 1.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0001.jpg

Figure 1

Extended structure of the cis-1,4-polyisoprene (PI) chain of pure PI terminated by H (_HPI_H), _ωPI_α6, and _ωPI_PO₄. End-isoprene residues of pure PI are shown by the cyan color. ω-Terminals (dimethyl allyl group and two trans-1,4-isoprene) are shown by blue color, α6 terminal is shown by orange color, and PO₄ is shown by red color. Backbone carbon and hydrogen atoms are shown in gray color and the oxygen and hydrogen atoms of –OH are shown in red and white colors, respectively.

Table 1

Details of the Three Melt Systems, _HPI_H, _ωPI_α6, and _ωPI_PO₄^a

system code	terminal-1	Cis-isoprene	terminal-2	N	M	NVT ns	simulated annealing ns	NPT ns	NVT production run ns
PI₀ (_HPI_H)	H-terminal	24	H-terminal	24	512	20	20	100	1000
PI_VI (_ωPI_α6)	ω	21	Hydroxy-terminated cis-1,4-isoprene (α6)	24	512	20	20	100	1000
PI_ph (_ωPI_PO₄)	ω	21	Phosphate-terminated cis-1,4-isoprene	24	512	20	20	100	5000

^aN is the number of monomers in each polymer chain and M is the number of chains present in each melt system. The full-atom MD simulations were performed in the sequence, i.e., NVT → simulated-annealing → NPT-EQ → NVT-production-run.

2. Methodology and Computational Discriptions

In this investigation, we explore three distinct types of cis-1,4-PI melt systems, each characterized by different combinations of ω-terminal and α′-terminal groups (details are provided in Table 1 and Figure Figure11). To ensure sufficient statistical sampling, we include 512 cis-1,4-PI chains in each melt system, with each chain comprising 24 monomers (M = 512 and N = 24). This approach aligns with our previous work, where we also considered a similar combination but solely focused on the six types of α-terminals.²²⁻²⁵ Experimental studies have revealed that NRs with low molecular weights exhibit enhanced mechanical and rheological properties, owing to a high density of terminal groups.⁶⁴⁻⁶⁶ These terminal groups readily associate, potentially interacting with nonrubber components such as proteins, lipids, and sugars, thereby facilitating the formation of the natural network structure.⁶⁶ In contrast, NRs with longer polymer chains exhibit poorer mechanical and physical properties due to the formation of a lower fraction of chain end group aggregates. Hence, we have chosen short PI chains in this study to gain molecular-level insights into terminal group aggregates. The molecular weight range of natural PI obtained from rubber trees spans from 100 kDa to 2000 kDa,^67,68 which is considerably large and too lengthy to be efficiently handled in numerical simulations. The primary objective of this study is to elucidate the difference in the formation of physical junction points among these terminals in the three types of polymer melt systems. Thus, we consider a cis-1,4-PI chain with a molecular weight of 1.6 kDa, which provides an adequate length for achieving the goals of this investigation.

2.1. All-Atom Molecular Dynamics Protocol

The initial configuration of each melt system is generated through a multiscale approach, as detailed in our previous publications.²²⁻²⁵ These initial structures (refer to Figure S1) are utilized as starting points for all-atom MD simulations, employing GROMACS version-2019.2.⁶⁹ The CHARMM general all-atom force field,^70,71 extensively validated for various polymer melt systems,⁶³ is employed for these simulations. The temperature and pressure are set to 360 K and 1 bar, respectively, following similar conditions from our previous work.^22,24,25 To prepare the systems, the initial structures are energy-minimized, and subsequently, NVT equilibration, NPT simulated annealing, and NPT equilibrations are performed (see Table 1 for the respective classical MD simulation times for each step). The equilibrated structure is then subjected to a production run spanning 1000–5000 ns. The trajectories generated from these simulations enable the computation of the radius of gyration, end-to-end distance, and radial distribution function between terminal groups, thereby elucidating the local structure of branching physical junction points between or among PI chain ends. To maintain the system temperature, we employ the Nosé–Hoover thermostat with a coupling constant of 1 ps. The system pressure is controlled at 1 bar using both the Berendsen barostat and Parrinello–Rahman barostat⁷² during equilibration and the production run, with time constants of 5 ps and a compressibility of 4.5 × 10^–5 bar^–1 (4.5 × 10^–10 Pa^–1), respectively. The Verlet cutoff scheme is employed to build the neighbor list, with a cutoff radius of 1.2 nm for all-atom MD simulations. The LINCS algorithm⁷³ is utilized to constrain HB lengths. A time step of 2 fs is employed for all-atom MD, while the particle mesh Ewald method⁷⁴ is employed for calculating electrostatic interactions.

3. Results

3.1. Equilibration of Polyisoprene Chains

In this section, we expound upon the outcomes derived from our all-atom MD simulations. The subject of our investigation revolves around the intricate assessment of static and dynamic properties, encompassing end-to-end distances, radii of gyration, as well as mean square internal distance and end-to-end vector autocorrelation functions, and Rouse mode analysis in relation to PI chains featured in the three distinct systems, as detailed in Table 1.

Through a systematic exploration, we scrutinized the equilibrium properties of PI chains, delving into key parameters such as the end-to-end distance (R_ee), radius of gyration (R_g), and mean square internal distances [left angle bracket] R²(n) [right angle bracket] . Capturing the temporal evolution of the end-to-end distance and radius of gyration via Figure S2 (provided in the Supporting Information), we observed remarkable stability over the course of the production run, signifying the establishment of equilibrium within the polymer chains across all PI melt systems. Moreover, we conducted a thorough evaluation, calculating the mean square end-to-end distance [left angle bracket] R_ee² [right angle bracket] and the mean square radius of gyration R_g², and the ensuing ratio R_ee²/R_g² is documented in Table 2. Our observations, encapsulated in Table 2, demonstrate the close adherence to a Gaussian behavior for the PI chains. The Kuhn length b characterizing a polymer is precisely defined as the quotient of the mean square end-to-end distance [left angle bracket] R_ee² [right angle bracket] and the maximum size when the polymer chain is fully extended R_max

Table 2

Mean-Square End-to-End Distance, R_ee², Mean Square Radius of Gyration, R_g², Kuhn Length, b, and Number of Kuhn Segments, N_RB, for all Melt Systems at 360 K and 1 bar

system	R_ee² [nm²]	R_g² [nm²]	R_ee²/R_g²	b [nm]	N_RB
_HPI_H	8.73 ± 0.25	1.43 ± 0.02	6.11	0.80	13.50
_ωPI_α6	9.00 ± 0.30	1.44 ± 0.03	6.23	0.81	13.59
_ωPI_PO₄	10.72 ± 0.15	1.61 ± 0.01	6.64	0.94	11.3

The variables b, representing the Kuhn length of the polymer chain, N_RB, denoting the number of Kuhn segments, and R_max, indicating the fully extended length of the polymer chain. Our investigation has involved the estimation of both the Kuhn length and the number of Kuhn segments for the polymer chain within each melt system, with detailed results presented in Table 2. The derived values for the Kuhn lengths exhibit close alignment with experimental data.⁷⁵

The investigation of chain conformations in each melt system was duly characterized through the utilization of mean square internal distances [left angle bracket] R²(n) [right angle bracket] , as elucidated by Figure S3 (contained in the Supporting Information). By averaging these quantities over all segments of size n = i – j along the chains, where i < j [set membership] [1, N] denotes monomer indices, we acquired profound insights into the overall system behavior. Notably, our findings underscore the sufficiency of the simulation duration for each melt system, as the chains exhibited extensive movement, spanning multiple times their individual dimensions.

The computation of the end-to-end vector autocorrelation function, denoted as C(t), was executed as an essential aspect of our investigation. This correlation function was assessed as a function of simulation time for each distinct melt system, employing the following well-established equation

Herein, R(0) and R(t) represent the end-to-end vectors at time t = 0 and time t, respectively, and the angle bracket signifies an ensemble average. The graphical depiction of C(t) for the three distinct systems (PI₀, PI_VI, and PI_ph) is presented in Figure Figure22. Additionally, we undertook an estimation of the average rotational relaxation time, achieved through a fitting procedure of the time correlation function to the simple exponential function as well as the stretched exponential function,⁷⁶ which is defined below

The rotational relaxation times, denoted as τ_rot, and the associated stretching exponents, represented as β for the three distinct systems, are illustrated in Figure Figure22b. In the case of the simple exponential function, β = 1 while for the stretched exponential function, β < 1. As seen from Figure Figure22a, the data for _HPI_H and _ωPI_α6 are well fitted by the simple exponential function, but for _ωPI_PO₄, the stretched exponential function gave a better fit. Intriguingly, the observed order of relaxation times exhibits a distinct hierarchy: _HPI_H < _ωPI_α6 [double less-than sign] _ωPI_PO₄. From the analysis of Figure Figure22b, we can see that the rotational relaxation times computed by using stretched exponential functions are always, respectively, smaller than the ones computed by using simple exponential decay functions. This compelling finding emphasizes that the presence of hydroxy and phosphate terminal groups significantly retards the dynamics of PI chains when compared to hydrogen-terminated PI chains. Especially the phosphate terminal groups alter the fashion of relaxation dynamics to the stretched exponential decay (β = 0.46) from a simple exponential decay one. In addition, in the _ωPI_PO₄ melt, the substantially long relaxation time τ_rot, (i.e., ten times or more, longer than the other two systems) is observed.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0002.jpg

Figure 2

Rotational correlation functions as a function of time for each melt system. In panel (a), the dashed lines represent fits to the stretched exponential function, An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m032.gif , on the other hand, the solid circle stands for a fit to a simple exponential function, C(t) = exp[−(t/τ_rot)]. Panel (b) showcases the corresponding rotational relaxation times, denoted as τ_rot. Computation of τ_rot is achieved by fitting the end-to-end vector autocorrelation functions, C(t), to a simple exponential function and stretched exponential function.

3.2. Rouse Mode Analysis

We conduct Rouse mode analysis on the all-atom chains, which are represented as coarse-grained chains comprising 12 Rouse beads based on the data, as shown in Table 2. In accordance with the Rouse model,⁷⁷ the normal coordinate is evaluated as X_p for free chain ends is defined as

Here, N_RB denotes the count of Rouse beads within a chain, while r_j denotes the position of the j-th Rouse bead.⁷⁸⁻⁸¹ In the context of a Rouse chain with free ends, the theoretical expectation dictates the relationship τ_p/τ_R = 1/p². Here, τ_R(≡τ₁) is the longest chain relaxation time.

The Rouse mode analysis for a polymer chain with fully immobilized chain ends can be conducted through the utilization of the subsequent normal coordinate⁸²

Under the condition of fixed chain ends, the anticipated theoretical relationship is expressed as τ_p/τ_R = 1/(p – 1/2)².

We computed the time autocorrelation function [left angle bracket] X_p(t)·X_p(0) [right angle bracket] /X_p(0)·X_p(0) for various Rouse modes, namely, p = 1, 2, 3, 4, 5, and 6, derived from the current all-atom MD simulations. The corresponding time autocorrelation functions are illustrated in Figure S4 of the Supporting Information. In Table 3, the longest chain relaxation time τ₁ corresponding to the first Rouse mode (p = 1), is depicted for every melt system. In Figure Figure22b, τ_rot., obtained through fitting a simple exponential decay function to the rotational correlation function C(t) for each melt system, is also presented. A consistent pattern is noted between the chain relaxation time of the first Rouse mode τ₁ (see Table 3) and the rotational relaxation time τ_rot. (Figure Figure22b).

Table 3

Chain Relaxation Time τ_p for First Rouse Mode (p = 1)^a

systems	type of alpha terminal	τ_R [ns]	τ_R (system)/τ_R (_HPI_H)
_HPI_H		48.88 ± 0.01	1.00
_ωPI_α6	hydroxy	57.07 ± 0.01	1.17
_ωPI_PO₄	phosphate	750.05 ± 0.25	15.34

^aThe fitting of the time autocorrelation function [left angle bracket]

X_p(t)·X_p(0) [right angle bracket]

X_p(0)·X_p(0) [right angle bracket]

with a simple exponential function for extracting the chain relaxation time τ_p.

Consistent with the predictions of the Rouse theory,⁷⁷ the envisaged scaling behavior of relaxation times τ_p corresponds to each normal mode p is expounded. The theoretical framework asserts that in the absence of constrained chain-ends, the scaling relationship is expressed as τ_p/τ_R = 1/p², while under fixed chain-ends, the relationship transforms to τ_p/τ_R = 1/(p – 1/2)².⁸²Figure Figure33 visually illustrates the dependence of the τ_p on the normal mode p. The determination of relaxation times τ_p for polymer chains in each melt system involves fitting the time autocorrelation function [left angle bracket] X_p(t)·X_p(0) [right angle bracket] /X_p(0)·X_p(0) with a simple exponential function, namely, exp(−t/τ_p). If Rouse scaling is valid, all curves presented in Figure Figure33 are anticipated to converge onto a unified trajectory, reflecting τ_p/τ_R = 1/p² for unconfined chain-ends and τ_p/τ_R = 1/(p – 1/2)² for fixed chain-ends. In the case of PI₀ (_HPI_H), where PI chains possess free chain-ends, adherence to Rouse scaling is evident (indicated by the gray color symbol ● in Figure Figure33 for Rouse modes p = 1, 2, 3, 4, 5, and 6). Likewise, for the _ωPI_α6 melt system, the PI chains also adhere to free-end Rouse scaling, suggesting that any deviation in higher modes is not triggered by intermittent associations among chain ends. In Table 3, we present the chain relaxation time τ_R for each molten system under consideration.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0003.jpg

Figure 3

Plot of τ_p vs p. Figure illustrates the relaxation times corresponding to distinct normal modes p of the polymer chains, elucidating the dynamic characteristics of the melt system. In the inset of the plot, we have shown the log–log plot of τ_p-vs-p.

The observations from Figure Figure33 and Table 3 collectively suggest that the system _ωPI_α6, characterized by end-group associations, maintains adherence to the scaling principles of the Rouse model with free ends. However, it manifests a prolonged relaxation time scale due to the intermittently constrained dynamics of chain ends. Conversely, the _ωPI_PO₄ melt system conforms similar behavior to 1/(p – 1/2)². In addition, from Table 3, we can see that the relaxation time can be substantially enhanced by PO₄ groups. Based on the findings of the current analyses, the intermittent interaction among α6-terminals does not alter the Rouse dynamics behavior; however, it does modulate the time scale of Rouse relaxation. Conversely, the association among PO₄ groups results in deviation from the behavior of a Rouse chain with free chain-ends, particularly for first and second Rouse modes (p = 1 and 2), and significantly augments the relaxation time of _ωPI_PO₄ chains.

3.3. Structure Aggregates of Hydroxylated- and Phosphorylated cis-1,4-Polyisoprene Chains

In this section, our focus centers on the examination of associations between terminal groups, specifically ISO–ISO, ω–ω, α6−α6, and PO₄–PO₄, in the PI₀, PI_VI, and PI_ph melt systems. Figure Figure44 presents snapshots of the PI₀, PI_VI, and PI_ph melt systems after the production run. Remarkably, we observed distinct characteristics in the aggregate states of α6−α6 and PO₄–PO₄, clearly evident in the snapshots, while such an aggregate state is notably absent for ISO–ISO.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0004.jpg

Figure 4

Snapshots of the (a) PI₀, (b) PI_VI, and (c) PI_ph melt systems captured after the final production-run of all-atom MD simulations. ω terminals are represented by the dimethyl allyl group (shown in blue), along with two trans-1,4-isoprene moieties. α6 and PO₄ terminals are depicted in orange and red colors, respectively, while the cis-1,4-isoprene units are represented in cyan. The backbone carbon and hydrogen atoms of each primary PI chain are illustrated in gray color, providing a comprehensive view of the molecular arrangement within the systems.

In order to gain insight into the local structural organization of terminal groups surrounding a specific terminal group, we have conducted a thorough investigation using radial distribution functions (RDFs). These RDFs, denoted as g_αβ(r), are estimated by assessing the distribution of distances between the centers of mass of various cis-1,4-PI chains’ terminal groups. The RDFs are mathematically described by the following equation

The expression for An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m007.gif represents the average density of type β particles in the vicinity of type α particles at a given distance r. This density is calculated as follows

In eq 7, An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m009.gif represents the local density of type β particles surrounding the i-th type α particle at a specific distance r. This local density is defined as follows

In eq 8, the primed summation notation introduces a distinction based on the relationship between α and β. Specifically, when α and β represent different particle types, the summation, denoted as ∑′, spans from j_β = 1 to N_β. Conversely, when α and β refer to the same particle type, the summation, represented as ∑′, excludes the contribution from the i_α-th particle and extends from j_β = 1 to N_β, with the exclusion of j_β ≠ i_α. The denominator on the right-hand side of eq 6, An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m011.gif , characterizes the average density of particle type β within a sphere of radius r_c, which is centered around the α particle. It is mathematically defined as follows

where [left angle bracket] ρ_β(r_c) [right angle bracket] _iα designates the density of particle type β confined within a sphere of radius r_c, centered around the i-th type α particle. This quantity is defined as follows

where V denotes the volume of the sphere, and the value of r_c is set to half of the box length. N_α represents the total number of α particles. Additionally, the Kronecker delta, δ_αβ, takes the values of 1 if α = β, and 0 if α ≠ β. To avoid double counting in the case α = β, the factor (δ_αβ + 1) is introduced in the denominator. In this context, our focus lies in the investigation of the solvation structure of terminal groups around a specific terminal group. This exploration offers valuable insights into the local molecular interactions and spatial arrangement of the terminal groups in the system under study. In the course of this investigation, a comprehensive analysis has been carried out to examine the RDFs pertaining to several key molecular components. Specifically, we have scrutinized the RDFs involving the center of mass of the dimethyl allyl ([DMA]) group, the center of mass of the H-terminated isoprene group ([ISO]), and the center of mass of the phosphate group (PO₄). Additionally, we have investigated the RDFs associated with the hydrogen atoms of the hydroxy groups within hydroxy-terminated cis-1,4-PI chains ([H]_OH) and the oxygen atoms of these same hydroxy groups ([O]_OH). These RDFs are denoted as [ISO]–[ISO], [O]_OH–[O]_OH, [O]_OH–[H]_OH [PO₄]–[DMA], [O]_OH–[DMA], and [DMA]–[DMA] within _HPI_H, _ωPI_α6, and _ωPI_PO₄ melt systems. For clarity, it is important to note that the notation [A]–[B] is employed to symbolically represent the arrangement of chemical functional group B around functional group A throughout this study. These RDFs are presented in Figure Figure55. In PI_0,VI,ph melt systems, the association between DMA groups exhibits a slightly higher degree compared to that of [ISO]–[ISO]. This observation aligns with our previous research findings.^22,23,25 In the case of PI_VI, the position of the first peak in the [O]-[H] RDF in _ωPI_α6 occurs at 0.2 nm, confirming the formation of HBs between the α6 and α6 terminal groups, which facilitates the formation of branch points between hydroxy-terminated PI chains. The intensity of this peak is significantly larger compared to the peaks in [ISO]–[ISO], [O]–[DMA], and [DMA]–[DMA] RDFs. This indicates a strong association of α6 terminal groups in PI chains compared to the ω terminals. In the case of _ωPI_PO₄, the [PO₄]–[PO₄] RDF peak intensity is notably higher than that of [PO₄]–[DMA] and [DMA]–[DMA] RDFs, implying the formation of a strong association between PO₄–PO₄ terminal groups. The number of large-sized polar aggregates between α6 and α6 and PO₄–PO₄ terminals, respectively, in _ωPI_α6 and _ωPI_PO₄ melt systems is significantly higher compared to [ISO]–[ISO] in _HPI_H. These findings shed light on the solvation structure and molecular associations of terminal groups in the PI systems under study.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0005.jpg

Figure 5

Radial distribution functions (RDFs) depicting the associations between terminal groups in the three types of cis-1,4-PI systems. For reference, the RDF of ISO around ISO in the pure PI system (PI₀) is plotted in all graphs, where ″[ISO]′′ refers to the center of mass of isoprene residues of terminals. Similarly, ′′[DMA]′′ denotes the center of mass of dimethyl allyl (dimethyl allyl residue of ω terminal) and ″[PO₄]″ represents the center of mass of the phosphate group. Moreover, [O] and [H] represent oxygen and hydrogen atoms of the hydroxy group in α6-terminals, respectively. In the RDFs, α and β are denoted by the symbols [PO₄], [ISO], [DMA], [O], and [H], and the corresponding associations α–β are color-coded as follows: gray for [ISO]–[ISO], blue for [DMA]–[DMA], red for [PO₄]–[PO₄], dark-green for [PO₄]–[DMA], black for [O]–[H], orange for [O]–[O], and dark-green for [O]–[DMA].

The peak intensity of [PO₄]–[PO₄] RDF is notably higher compared to [α6]–[α6] and [ISO]–[ISO] RDFs (see Figure S5). This observation indicates that the [PO₄] terminal groups of PI chains exhibit stronger associations than the hydroxy-terminated PI chains. Notably, we observed narrower RDF peaks between the phosphate terminal groups of phosphate-terminated PI chains, while the first RDF peaks of α6 terminal groups of hydroxy-terminated PI chains are broader. Consequently, the first coordination shell of PO₄ terminals appears to be more ordered than that of α6 terminals. Additionally, the coordination shell of α6 around α6 is more structured compared to that of ISO around ISO. The analysis of running coordination numbers (RCNs) (Figure S6) further supports these findings. Specifically, in the first coordination shell, the number of PO₄ terminal groups around other PO₄ terminals is significantly larger compared to the α6 terminals around other α6 terminals and ISO around ISO terminals. Notably, the value of RCNs in the first coordination shell is highest for PO₄ terminal groups, providing crucial insights into the distinct intermolecular associations and structural order of terminal groups in the PI systems under study.

To investigate the formation of clusters by terminal groups, we employed the single linkage algorithm for association analysis as described in refs (83–85) In the subsequent analysis, we define an “encountering state” of an end-group of a chain with other end-groups that belong to different PI chains. Specifically, when the distance between two end-groups from different chains falls below a threshold distance r_th at a specific time step, we identify these end-groups as being in an encountering state. Moreover, if either of the two end-groups is in the encountering state with a different end-group from the other end, we consider the third end-group as being in the encountering state with the first two end-groups. By iteratively applying this procedure to newly added end-groups and stopping when there are no further end-groups to be added, we can determine the size s of the encountering end-groups by counting the number of end-groups involved. Consequently, we evaluate the number n_s of encountering states with size s, thus providing valuable insights into the clustering behavior of terminal groups in the PI system under investigation.

The classification of end-groups in an encountering state is based on whether they contain only a single type of end-group, denoted as α type, in which case the notation n_s^αα is employed, or whether they involve two different types of ends, α and β, in which case the notation n_s^αβ is used. To evaluate the number of encountering end-groups of size s, we perform this analysis for K instances during the last 100 ns of the final production run trajectories. Here, K represents the number of frames present in the 100 ns all-atom simulation trajectories, i.e., n_s^(αβ)(k), where k = 1, 2, ..., K = 50,000. Using n_s^(αβ)(k), we define the encountering-event-fraction f_enc^(αβ)(s) for a given cluster size s as follows

where the number of chains in the system is denoted as M. In Figure Figure66, we present f_enc^(αβ)(s) as a function of the encountering-event size s for each melt system. Specifically, f_enc^(αβ)(s = 2) represents the probability of a terminal group encountering another terminal group, effectively indicating the likelihood of terminal group associations in the form of dimers. If f_enc^(αβ)(s = 2) = 1, it implies that every terminal group from any chain must be associated with a terminal group from another chain, resulting in 100% of terminal groups forming dimers. However, as the size of the encountering-event increases, the encountering-event-fraction of terminal groups decreases. This suggests that the probability of terminal group associations decreases with the increasing size of the encounter, highlighting the diversity in the sizes and structures of the associations in the PI melt systems. In the case of _HPI_H, our observations reveal that the encountering-end-groups predominantly consist of size two, with a few occurrences of size three, and very few instances of size four. However, for _ωPI_α6, the situation is different. Here, in addition to sizes two and three, encountering-events of larger sizes beyond three is also present due to the strong attractive interactions between α6 and α6 terminal groups, as corroborated by the findings presented in Figure Figure55. The clusters of size s are evident as (OH)_pp(DMA)_qq or An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m015.gif , where pp = 0, 1, 2,···and qq = 0, 1, 2,···with s = pp + qq. In terms of the encountering-event-fractions of size two, the encounter of [(OH)₂] is slightly more frequent than [(OH)₁(DMA)₁] and [(DMA)₂]. Moving to size three, [(OH)₃] also dominates over [(DMA)₃], [(OH)₂(DMA)₁], and [(OH)₁(DMA)₂]. In the case of size four, [(OH)₄] significantly participates, while the contribution of [(OH)₁(DMA)₃] is minor (the presence of other possibilities like [(OH)₂(DMA)₂ is not observed). For size five, only [(OH)₅] is observed, and for size six, only [(OH)₆] exists. The encountering-event-fractions in the _ωPI_α6 melt systems are more pronounced compared to the _HPI_H system. In the _ωPI_PO₄ system, the encountering-event-fractions of size two are dominated by [(DMA)₂], with the contribution of An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m016.gif being very small, while is almost absent. Regarding size three, [(DMA)₃] is more widespread compared to and while the contribution of is negligible. For the encountering-event-fractions of size four, [(DMA)₄] predominates over An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m021.gif (the presence of other possibilities like is not observed). For encountering-event-fractions of size larger than four (4 [less-than-or-eq, slant] s 28), exclusively exists. The observation indicates that the end-groups encounter each other with a certain fraction; however, their stability as branch points or stable clusters remains undetermined. To address this, potentials of mean force (PMFs) have emerged as a widely used approach for assessing the stability of clusters in various research studies.^{25,62,86−101} In this study, we compute the PMFs between the terminal groups using the following equation

where T is the temperature of the system, k_B is the Boltzmann constant in kJ mol^–1/K, and g(r) represents the radial distribution function between the terminal groups. This equation allows us to derive the PMFs, which provide essential insights into the thermodynamic stability and interactions of the terminal groups in the system. We present the PMF between the terminal groups as a function of distance in Figure Figure77. Each PMF profile exhibits a minimum for contact pairs, which we denote as the contact minimum (CM). In PI₀, we observe a shallow minimum for contact pairs between the two isoprene terminals, [ISO]-[ISO], at a distance of 0.6 nm. Similarly, for PI₀ and PI_VI melt systems, we perceive a shallow minimum for contact pairs between DMA- and DMA-terminals at 0.6 nm. The interaction free energy of these contact pairs is on the order of thermal energy, indicating that they cannot be considered stable associations. In contrast, for the _ωPI_α6 system, we observe sharp minima for contact pairs between [OH]–[OH] of α6 terminals, and the stability of this contact pair is substantially higher than the thermal energy. The stability of α6–α6 association is significantly higher than that of [ISO]–[ISO] and [DMA]–[DMA] terminal associations. In the case of the _ωPI_PO₄ system, we notice the formation of a stable contact pair between [PO₄]–[PO₄] at a distance of 0.44 nm and between [DMA]–[DMA] at 0.6 nm. The [PO₄]–[PO₄] contact pair is more stable than the [DMA]–[DMA] contact pair. For the formation of a stable cluster, the absolute value of the interaction energy |W| at the CM must be significantly greater than the thermal energy. Therefore, we define the criteria for cluster formation as |W| > 2k_BT at the CM. These PMFs provide valuable information about the stability of terminal group associations and are crucial for understanding the behavior of the terminal groups in the PI systems studied. The dissociation barrier height for the CM of end groups is defined as the free energy difference between the CM and the first maximum of the PMF.¹⁰²Figure Figure77 clearly shows that the values of dissociation barriers are 10.8 kJ/mol (3.6 k_BT) for the [O]-[O] CM and 17.9 kJ/mol (5.9 k_BT) for the [PO₄]–[PO₄] CM in the _ωPI_α6 and _ωPI_PO₄ melt systems, respectively. These dissociation barrier heights indicate the energy required to break the stable contact pairs, and they offer insights into the stability of terminal group associations in the respective PI systems.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0006.jpg

Figure 6

Plot illustrates the fraction of end-group encounters with a cluster size s for the PI_0,VI,ph systems. Cut threshold length, r_th, is set to 0.6 nm. It is important to note that the sum of the fractions for all cluster sizes, ∑_s = 1^∞f_enc^(αβ) (s = 1), is equal to 1 for each system, ensuring a comprehensive account of all cluster sizes. However, the fractions corresponding to s = 1 are not explicitly shown in the plot.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0007.jpg

Figure 7

PMF W(r) is calculated by using the equation W(r) = −k_BT log(g(r)), where g(r) is the radial distribution function between the terminal groups. Error bars in each PMFs profile is less than 0.2 k_BT.

By utilizing the PMF profiles between the terminal groups, we define a stable association between two terminal groups as follows: when the PMF (W) between two terminal groups satisfies the inequality (|W(r_CM)| > 2k_BT), where r_CM is the distance at the CM, the two terminal groups are considered to be in a state of stable association and can form a “cluster”. To count the number of terminal groups that join a cluster, we apply two criteria: (i) the distance criterion mentioned in eq 11 and (ii) the criterion for stable association (|W(r_CM)| > 2k_BT). The obtained number of clusters with size s is then expressed as n_s^(αβ)(k;|W(r_CM)|>2k_BT), where α and β represent OH or PO₄ or DMA. The meaning of the superscripts α and β is the same as that used for the number of encountering end groups. For example, An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m025.gif is the number of size s clusters that contain PO₄ and DMA residues of the terminals, and n_s^(OH–OH) (n_s^(DMA–DMA) or ) is the number of size s clusters composed of only OH (PO₄ or DMA) terminal groups. By using these numbers, we determine the cluster-formation-fraction f_cluster^(αβ)(s) with the following equation

The cluster-formation-fraction, f_cluster^(αβ)(s), is presented in Figure Figure88. In the case of _HPI_H, the formation of clusters by the terminal groups is entirely absent, while the encountering-events of end-groups are clearly evident, as illustrated in Figure Figure66. Conversely, for _ωPI_α6, clusters of sizes two, three, four, five, and six are observed, with no clusters of size s ≥ 7 being detected. In the context of _ωPI_PO₄, clusters of size four, and larger than four are identified. Notably, the associations involving [DMA]–[DMA] and [OH]–[DMA] terminal groups are entirely excluded from cluster formation, aligning with a similar trend reported in our previous work.²⁵ Furthermore, in the case of _ωPI_PO₄, cluster-formation-fractions are observed for 4 ≤ s ≤ 28, with the associations of [DMA]–[DMA] and [PO₄]–[DMA] terminal groups being fully precluded from cluster formation.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0008.jpg

Figure 8

Cluster-formation-fraction as a function of cluster size s in PI_0,VI,ph systems. Additional criterion utilizing the magnitude of PMF between terminal groups to determine stable cluster formation. Note that ∑_s = 1^∞f_cluster^(αβ)(s) = 1 for each system, and the fractions for s = 1 are not shown.

3.4. Survival Time of Terminal Groups in the Coordination Shell

The local rotational and translational mobility of terminal groups is influenced by the duration of time they spend within a specific coordination shell. The survival probability P(τ) quantifies the fraction of terminal groups that remain within the coordination shell around a given terminal group at time t and continue to be present in the same shell at a later time t + τ. This probability is determined by assessing the likelihood of finding a group of particles remaining in a spherical region with a radius of r_thrd, centered at the geometrical center of the selected group.^103,104 The survival probability is calculated using the following equation

where N(t) is the number of terminals in a spherical region with a radius (r_thrd) centered at the geometrical center of the selected group assigned to a cluster at time t, N(t, t + τ) is the number of terminals that continue to stay in the same spherical region from time t to a time t + τ, M_tot is the total number of time steps contributing to P(τ), and τ is the time between the analyzed configurations. The radius r_thrd of the selected region for the calculation of P(τ) for ISO around ISO, DMA around DMA, OH around OH, and PO₄ around PO₄ is chosen based on the RDFs between [ISO]–[ISO], [DMA]–[DMA], [O]–[O], and [PO₄]–[PO₄] (Figure Figure55). The radius of the first coordination shells of ISO around ISO, DMA around DMA, and O around O is used for the selected regions. MDAnalysis¹⁰⁵ is employed for the calculation of the survival probability. From the investigation of Figure Figure99, the order of decay of survival probability is _HPI_H [double less-than sign] _ωPI_α6 _ωPI_PO₄. The survival probabilities slowly decay in phosphate-terminated PI chains compared to hydroxy-terminated PI. The decay times of the survival probability curves are obtained by fitting them into the Kohlrausch–Williams–Watts (KWW) stretched exponential function, and the estimated values are 6.05 ps _HPI_H, 56.22 ps (_ωPI_α6), and 462 ps (_ωPI_PO₄). The fitted data is shown in Figure S7 of the Supporting Information. These values indicate a much slower exchange rate of terminal groups within the first coordination shell in both phosphorylated and hydroxylated PI chains compared to pure PI. Furthermore, the exchange rate of PO₄ terminal groups within the first coordination shell is much slower compared to hydroxy terminal groups. The decay rates and stretching exponents, along with their error bars, are presented in Table 4.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0009.jpg

Figure 9

Plot shows the survival probabilities P(t) for terminal groups within the first coordination shell of each terminal group as functions of time duration of t. Survival probabilities were calculated using the last 10 ns of the 1000 ns independent MD simulations. Error bars were calculated using the block-average method. Each survival probability was computed by averaging over the available time windows at given intervals ranging from 1 to 5000 ps.

Table 4

Survival Probability P(t) Is Analyzed Using the Kohlrausch–Williams–Watts Stretched Exponential Function, Enabling the Computation of the Average Relaxation Time τ and the Stretching Exponent β

systems	average relaxation time [ps]	stretching exponents β
_HPI_H	6.05 ± 0.05	0.668 ± 0.001
_ωPI_α6	56.22 ± 0.25	0.665 ± 0.002
_ωPI_PO₄	462 ± 4	0.143 ± 0.001

The relaxation time of phosphate terminal groups interacting with other phosphate terminal groups is considerably longer compared to the relaxation time of hydroxy terminal groups interacting with other hydroxy terminal groups.

3.5. Self-Diffusion Behavior

We analyze the diffusion behavior of individual PI chains in each melt system. The mean square displacement (MSD) of the center of mass of a single PI chain is computed using the last 100 ns of the production-run trajectories

where ΔR_cm(t) ≡ R_cm(t) – R_cm(0) and [left angle bracket] ΔR_cm²(t) [right angle bracket] represent the MSD of the center of mass of a PI chain at time t. The resulting MSD as a function of simulation time is presented in Figure Figure1010. As seen in Figure Figure1010, the MSD does not reach the linear regime yet, but we can see that the self-diffusion of PI with the phosphate terminal is suppressed and the diffusion exponent α is substantially smaller than that of _HPI_H and _ωPI_α6. We observed that the center of mass motion of polymer chains is subdiffusive for all three melt systems, i.e., [left angle bracket] ΔR_cm²(t) [right angle bracket] grows as t^0.70, t^0.69, and t^0.49 for _HPI_H, _ωPI_α6, and _ωPI_PO₄, respectively. A similar trend was reported in previous experimental and theoretical studies.^106−110

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0010.jpg

Figure 10

Mean square displacements for the centers of mass of polymer chains in each melt system. Inset stands for the log–log plot of the main graph.

We observed that the self-diffusion follows the order: _HPI_H > _ωPI_α6 [dbl greater-than sign] _ωPI_PO₄. Thus, the single-chain mobility of the pure-PI system is faster compared to _ωPI_α6 and _ωPI_PO₄ melt systems. The presence of HBs between α6−α6 and PO₄–PO₄ terminals appears to retard the mobility of hydroxy-terminated and phosphate-terminated PI chains. In particular, we found that the diffusion of _ωPI_PO₄ molecular chains is significantly slower than that of _ωPI_α6 molecular chains.

3.6. Characterization of Polar Physical Junction Points: Insights into Structural Terminal Groups

The investigation into the dynamics of HBs within the polymer melt systems involved the computation of the widely recognized HB time correlation function C_HB(t) for pairs of hydrogen-bonded polymer chains denoted as i, j. This correlation function, essential for understanding HB dynamics, follows the definition established in previous studies^111−114

Here, the dynamic nature of HBs in polymer chains is examined through the variable h_ij(t), which quantifies the fulfillment of geometric HB criteria between pairs of polymer chains i and j at time t. The summation covers all potential pairs exhibiting hydrogen-bonding interactions. The angular brackets represent an average over multiple starting times within the trajectory. In the continuous HB correlation framework C_HB^c(t), it is crucial to highlight that a HB, once disrupted, maintains its broken classification, irrespective of any subsequent reformation. This definition offers insights into the average duration a pair remains bonded, providing the average HB lifetime, denoted as the continuous lifetime τ_HB^c. In contrast, intermittent HB correlation C_HB^I(t) explores the persistence probability of a HB formed at t = 0, despite multiple breakings and reformations within the time interval [0, t]. The corresponding lifetime is referred to as the intermittent lifetime or HB relaxation time τ_HB^I. We computed both continuous C_HB^c(t) and intermittent C_HB^I(t) HB correlation functions for α6−α6 and PO₄–PO₄ utilizing MDAnalysis.^105,115 We conducted an assessment of the continuous C_HB^c(t) and intermittent C_HB^I(t) HB correlation functions for specific pairs, including α6−α6 and PO₄–PO₄, as illustrated in Figure Figure1111. Furthermore, the determination of the HB relaxation time τ_HB was achieved through fitting the HB autocorrelation functions C_HB^c(t) and C_HB^I(t) using the KWW stretched exponential function. The resultant values of τ_HB and the stretching exponent β_HB are documented in Table 5. For the intermittent HB correlation function, we employed the KWW stretched exponential function, specifically expressed as An external file that holds a picture, illustration, etc.
Object name is lg4c00019_m031.gif . The estimated values of τ_HB^I, stretching exponent β_HB^I, and P_HB^I are presented in Table 5, offering a detailed characterization of the temporal behavior and persistence of HBs within our studied systems.

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0011.jpg

Figure 11

Hydrogen bond correlation functions, the continuous C_HB^c(t) (A), and intermittent C_HB^I(t) (B), for pairs α6−α6 (orange) and PO₄–PO₄ (red).

Table 5

HB Relaxation Time τ_HB and Stretching Exponent β_HB for _ωPI_α6 and _ωPI_α6 Melt Systems

	continuous		intermittent
system	τ_HB^c [ps]	β_HB^c	τ_HB^I [ps]	β_HB^I	P_HB^I
_ωPI_α6	11.50 ± 0.01	0.798 ± 0.001	32.79 ± 0.18	0.2697 ± 0.0007	0.0170 ± 0.0002
_ωPI_PO₄	522.40 ± 0.66	0.610 ± 0.001	371,193 ± 10,720	0.513 ± 0.005	0.0738 ± 0.0007

Examination of the continuous and intermittent HB lifetimes, as elucidated in Figure Figure1111 and summarized in Table 5, underscores a pronounced distinction between the interactions involving PO₄–PO₄ and α6−α6 pairs. The notably protracted continuous lifetime exhibited in PO₄–PO₄ interactions signifies a considerably more robust correlation among hydrogen-bonded _ωPI_PO₄ PI chains in comparison to their _ωPI_α6 counterparts. This observation suggests a higher degree of order and structural organization within the polar physical junctions containing phosphate groups as opposed to those involving OH terminal groups.

4. Discussion and Concluding Remarks

We have conducted all-atom MD simulations to investigate three distinct types of cis-1,4-PI melt systems featuring various combinations of ω and α′ terminal groups. Our analysis encompasses both static and dynamic characteristics of PI chains, including end-to-end distance R_ee, radius of gyration R_g, mean square internal distances, end-to-end vector autocorrelation function C(t), average relaxation time τ_rot, Rouse mode analysis, survival probability P(t) of terminal groups, HB correlation functions of terminal groups, and self-diffusion behavior of PI chains. The examination of the end-to-end vector autocorrelation function depicted in Figure Figure22a has unveiled a discernible and intriguing behavior: a sluggish relaxation pattern is evident in PI chains featuring PO₄ terminal groups. This tardy mobility is further substantiated by the reduced translational motion of PI chains in the context of the PO₄-terminated PI melt system, as shown in Figure Figure1010. This behavior strongly implies the likelihood of clustering phenomena between phosphate-terminated PI chains. The Rouse mode analyses (Figure Figure33) reveal that the intermittent interaction among α6-terminals does not induce alterations in Rouse dynamics behavior; nevertheless, it does exert modulation on the time scale of Rouse relaxation. However, the _ωPI_PO₄ melt system showed similar behavior to 1/(p – 1/2)², concurrently resulting in a substantial increase in the relaxation time of _ωPI_PO₄ chains. We noted a congruent trend in the observed patterns between the chain relaxation time of the first Rouse mode τ₁ (Table 3) and the rotational relaxation time τ_rot (Table of Figure Figure22). The RDFs between various terminal groups, namely, [ISO]–[ISO], [DMA]–[DMA], [O]_OH–[O]_OH, [α6]–[α6], and [PO₄]–[PO₄] (as illustrated in Figures Figures55 and S5), offer additional valuable insights into these interactions. For [ISO]–[ISO] and [DMA]–[DMA], the low intensity and broad first RDF peaks indicate a relatively weak association between these groups. Conversely, the presence of a sharper and relatively higher first RDF peak for [DMA]–[DMA] than that of [ISO]–[ISO] suggests a more stable and structured association among ω terminals, especially when compared to the isoprene terminals in the pure PI system PI₀. In the context of _ωPI_PO₄, the most remarkable observation is the significantly higher intensity of the first peak in the [PO₄]–[PO₄] RDFs, in contrast to [ISO]–[ISO] in pure PI, [O]_OH–[O]_OH in _ωPI_α6, and ω–ω in PI_0,VI,ph systems. This heightened intensity underscores a much stronger association between PO₄ and PO₄ terminals in the specific context of the _ωPI_PO₄ system. This significant finding is further substantiated by the slower decay of the survival probability for PO₄ around PO₄ terminals, as demonstrated in Figure Figure99. It is crucial to emphasize that this robust association between phosphate groups is pivotal in the formation of a naturally occurring network, a trend that aligns with prior experimental investigations.³³ This intriguing phenomenon sheds light on the intricate MD and intermolecular interactions in the PI system, especially in the presence of phosphate-terminated chains, which have profound implications for material properties and behavior.

The investigation into the clustering and dynamics of terminal groups within physical junction points was conducted through a comprehensive analysis of three distinct PI melt systems. Our focus centered on assessing the encountering-event-fraction f_enc^(αβ)(s) of terminal groups (Figure Figure66), the PMF between these terminal groups (Figure Figure77), and the cluster-formation fraction of terminal groups (Figure Figure88). The outcomes of these analyses provided valuable insights into the unique characteristics exhibited by each melt system. The _ωPI_PO₄ melt system demonstrated a prevalence of encountering events of larger sizes, specifically over the range of (4 [less-than-or-eq, slant] s 28), whereas the _HPI_H and _ωPI_α6 systems exhibited encountering events of smaller sizes (2 ≤ s ≤ 6). The structural analysis of the _ωPI_PO₄ system further revealed a sharp CM in the PMF W(r) at approximately 0.44 nm. Importantly, the depth of W(r_CM) at this CM significantly exceeded the thermal energy (2.99 kJ/mol at 360 K). This observation indicated that the association of PO₄ terminals was notably more stable than the associations of α6−α6 in PI_VI and those of ω–ω in PI_VI and PI_ph. This finding was further substantiated by the survival probabilities of terminal groups, reinforcing the enhanced stability of PO₄ associations. Our systematic analysis not only elucidates the distinctive clustering behavior and dynamics within the studied PI melt systems but also highlights the stabilizing effect of phosphate-terminated PI chains, contributing valuable insights to the understanding of molecular interactions in polymer systems.

Our investigation delved into a comprehensive analysis of the cluster-formation-fraction f_cluster^(αβ)(s), enriched by the incorporation of an additional criterion, (|W(r_CM)| > 2k_BT). In the _ωPI_PO₄ system, a noteworthy revelation unfolded as larger and remarkably stable clusters, spanning the range of four to twenty-eight, emerged. These clusters served as clear indicators of substantially expanded physical junction points within the system. This stands in marked contrast to the _ωPI_α6 system, where hydroxy terminal groups orchestrated the establishment of more compact physical junction points, varying in size from two to seven, as depicted in Figure Figure88. Visual representations capturing the essence of these physical junction points are presented in Figure Figure1212, complemented by Videos S1 and S2 for the _ωPI_α6 and _ωPI_PO₄ melt systems, respectively. An intriguing alignment with recent experimental findings was observed, underscoring the pivotal role played by phosphate terminal groups in orchestrating the formation of extensive polar aggregates, synonymous with physical junction points, in phosphate-terminated PI chains, as documented by Li et al.³³ In the context of the pure PI melt system, cluster formation exhibited complete suppression due to the feeble association between isoprene-terminal groups.²²⁻²⁵

An external file that holds a picture, illustration, etc.
Object name is lg4c00019_0012.jpg

Figure 12

Visualization of physical junction points formed by α6 and PO₄ terminals within PI_VI (a) and PI_ph (b) melt systems, respectively. Dimethyl allyl group of ω terminals is represented in blue, while α6 and PO₄ terminals are color-coded in orange and red, respectively. Depicted orange and red regions are selectively chosen around the atoms within the physical junction points, with a radius of 0.8 nm. Graphical representation employs the quicksurf drawing method with a radius scale of 1.5, a density isovalue of 9.0, and a grid spacing of 2.0. Backbone carbon and hydrogen atoms are illustrated in gray for clarity and context.

To recapitulate, our investigation probed the impact of ω and α′ terminals on both the association of end-terminal groups and the dynamic behavior of PI chains within three distinct melt systems. The presence of a high intensity first RDF peak at 0.2 nm between the [O]_H₂PO₄–[H_(OH)1]_H₂PO₄ and [O]_H₂PO₄–[H_(OH)2]_H₂PO₄ groups confirms the formation of a strong HB between the H₂PO₄ groups (Figure S8). The robust hydrogen-bond interactions observed in H₂PO₄–H₂PO₄ terminal pairs (Figure S8) facilitated the formation of substantial, stable clusters with sizes greater than 2. These polar aggregates, acting as dynamic cross-linking sites, heightened the cross-linking density, thereby amplifying the mechanical properties of the material. In contrast, the comparatively weaker associations between [ISO]–[ISO] and ω–ω terminal pairs precluded the formation of stable clusters, resulting in the absence of steadfast branching points. Consequently, no enduring junctions emerged between ISO–ISO in the case of _HPI_H and between ω–ω in the scenarios of _ωPI_α6 and _ωPI_PO₄. This detailed investigation illuminates the subtle interplay between terminal associations and cluster stability, offering valuable insights into the mechanical properties of the examined PI melt systems. We have validated the presence of highly stable polar aggregates, formed by robust hydrogen-bond interactions among the phosphate terminals of PI chains. These aggregates hold the potential to be a key factor contributing to the occurrence of SIC observed in NR.

In future research pursuits, we intend to expand our investigation through comprehensive all-atom MD simulations applied to hydroxylated- and phosphorylated-PI melt systems. Our approach will extend beyond traditional rubber constituents to include nonrubber elements such as proteins and phospholipids. This extended exploration seeks to elucidate the intricate interplay between nonrubber components and the formation of aggregates, whether polar or nonpolar, involving hydroxy- and phosphate-terminal groups of PI chains. The anticipated findings hold the potential to provide deeper insights into the complex molecular interactions governing the exceptional properties of NR, particularly its heightened toughness and the intriguing phenomenon of SIC.

Acknowledgments

This work was supported by JST and CREST grant no. JPMJCR2091, Japan. Our heartfelt appreciation extends to Professor Kenji URAYAMA from the Department of Material Chemistry, Graduate School of Engineering, Kyoto University, whose invaluable insights and intellectually stimulating discussions have notably enhanced the outcomes of this research endeavor. We thank Shoma Fujii for providing his code for mean square internal distance and Rouse mode analysis. Furthermore, we extend our gratitude to the Center for Computational Materials Science at the Institute for Materials Research, Tohoku University, for granting us access to the MASAMUNE-IMR platform through project nos. 202112-SCKXX-0017, 202212-SCKXX-0012, and 202312-SCKXX-0004. We also convey our sincere appreciation to the “Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures” and the “High-Performance Computing Infrastructure” in Japan (ProjectID nos. jh210017-MDH, jh220054, jh230061, and jh240063) for their invaluable provision of computational resources. The utilization of computational systems such as Wisteria/BDEC-01, hosted at the Information Technology Center, University of Tokyo, along with Octopus and SQUID located at the Cybermedia Center, Osaka University, has significantly contributed to the successful execution of this research.

Supporting Information Available

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acspolymersau.4c00019.

Initial structures of PI melt systems; end-to-end distance (R_ee); radius of gyration (R_g); mean square internal distances; time autocorrelation function of normal coordinate of different Rouse modes; RDFs; RCN; and fitted survival probabilities P(t) (PDF)
Formation of polar aggregates between α6 terminals in the _ωPI_α6 melt system (MP4)
Visual depiction of the formation of polar aggregates between PO₄ terminals in the _ωPI_PO₄ melt system (MP4)

Notes

This work was supported by JST and CREST (grant no. JPMJCR2091), Japan.

Notes

The authors declare no competing financial interest.

Supplementary Material

lg4c00019_si_001.pdf^{(14M, pdf)}

lg4c00019_si_002.mp4^{(52M, mp4)}

lg4c00019_si_003.mp4^{(51M, mp4)}

References

Oouchi M.; Ukawa J.; Ishii Y.; Maeda H. Structural Analysis of the Terminal Groups in Commercial Hevea Natural Rubber by 2D-NMR with DOSY Filters and Multiple-WET Methods Using Ultrahigh-Field NMR. Biomacromolecules 2019, 20, 1394–1400. 10.1021/acs.biomac.8b01771. [Abstract] [CrossRef] [Google Scholar]
Tarachiwin L.; Sakdapipanich J.; Ute K.; Kitayama T.; Bamba T.; Fukusaki E. I.; Kobayashi A.; Tanaka Y. Structural characterization of α-terminal group of natural rubber. 1. Decomposition of branch-points by lipase and phosphatase treatments. Biomacromolecules 2005, 6, 1851–1857. 10.1021/bm058003x. [Abstract] [CrossRef] [Google Scholar]
Tanaka Y. Structural characterization of natural polyisoprenes: Solve the mystery of natural rubber based on structural study. Rubber Chem. Technol. 2001, 74, 355–375. 10.5254/1.3547643. [CrossRef] [Google Scholar]
Cornish K. Similarities and differences in rubber biochemistry among plant species. Phytochemistry 2001, 57, 1123–1134. 10.1016/S0031-9422(01)00097-8. [Abstract] [CrossRef] [Google Scholar]
van Beilen J. B.; Poirier Y. Establishment of new crops for the production of natural rubber. Trends Biotechnol. 2007, 25, 522–529. 10.1016/j.tibtech.2007.08.009. [Abstract] [CrossRef] [Google Scholar]
Le Cam J. B.; Toussaint E. The mechanism of fatigue crack growth in rubbers under severe loading: The effect of stress-induced crystallization. Macromolecules 2010, 43, 4708–4714. 10.1021/ma100042n. [CrossRef] [Google Scholar]
Toki S.; Burger C.; Hsiao B. S.; Amnuaypornsri S.; Sakdapipanich J.; Tanaka Y. Multi-scaled microstructures in natural rubber characterized by synchrotron X-ray scattering and optical microscopy. J. Polym. Sci., Part B: Polym. Phys. 2008, 46, 2456–2464. 10.1002/polb.21578. [CrossRef] [Google Scholar]
Gent A. N.; Kawahara S.; Zhao J. Crystallization and Strength of Natural Rubber and Synthetic cis-1,4-Polyisoprene. Rubber Chem. Technol. 1998, 71, 668–678. 10.5254/1.3538496. [CrossRef] [Google Scholar]
Amnuaypornsri S.; Toki S.; Hsiao B. S.; Sakdapipanich J. The effects of endlinking network and entanglement to stress-strain relation and strain-induced crystallization of un-vulcanized and vulcanized natural rubber. Polymer 2012, 53, 3325–3330. 10.1016/j.polymer.2012.05.020. [CrossRef] [Google Scholar]
Tang M.; Zhang R.; Li S.; Zeng J.; Luo M.; Xu Y. X.; Huang G. Towards a Supertough Thermoplastic Polyisoprene Elastomer Based on a Biomimic Strategy. Angew. Chem., Int. Ed. 2018, 57, 15836–15840. 10.1002/anie.201809339. [Abstract] [CrossRef] [Google Scholar]
Payungwong N.; Inoue T.; Sakdapipanich J. A Model Study of the Influence of the Natural Rubber (NR)- Endogenous Gel Fraction on the Rheological Performance of NR Using Synthetic Polyisoprene Rubber (IR) Blends with Different Ratios of Gel. ACS Appl. Polym. Mater. 2022, 4, 7061–7069. 10.1021/acsapm.2c00979. [CrossRef] [Google Scholar]
Wu J.; Qu W.; Huang G.; Wang S.; Huang C.; Liu H. Super-Resolution Fluorescence Imaging of Spatial Organization of Proteins and Lipids in Natural Rubber. Biomacromolecules 2017, 18, 1705–1712. 10.1021/acs.biomac.6b01827. [Abstract] [CrossRef] [Google Scholar]
Toki S.; Che J.; Rong L. S.; Hsiao B.; Amnuaypornsri S.; Nimpaiboon A.; Sakdapipanich J. Entanglements and Networks to Strain-Induced Crystallization and Stress-Strain Relations in Natural Rubber and Synthetic Polyisoprene at Various Temperatures. Macromolecules 2013, 46, 5238–5248. 10.1021/ma400504k. [CrossRef] [Google Scholar]
Kawahara S.; Matsuura A.; Kakubo T.; Nishiyama N.; Tanaka Y. Crystallization Behavior of Natural Rubber, Part 3. Effect of Saturated Fatty Acid on Crystallization of cis-1,4-Polyisoprene. Nippon Gomu Kyokaishi 1999, 72, 413–417. 10.2324/gomu.72.413. [CrossRef] [Google Scholar]
Kawahara S.; Tanaka Y. Plasticization and crystallization of cis-1,4 polyisoprene mixed with methyl linolate. J. Polym. Sci., Part B: Polym. Phys. 1995, 33, 753–758. 10.1002/polb.1995.090330503. [CrossRef] [Google Scholar]
Chen Q.; Zhang Z.; Huang Y.; Zhao H.; Chen Z.; Gao K.; Yue T.; Zhang L.; Liu J. Structure-Mechanics Relation of Natural Rubber: Insights from Molecular Dynamics Simulations. ACS Appl. Polym. Mater. 2022, 4, 3575–3586. 10.1021/acsapm.2c00147. [CrossRef] [Google Scholar]
Xu L.; Huang C.; Luo M.; Qu W.; Liu H.; Gu Z.; Jing L.; Huang G.; Zheng J. A rheological study on non-rubber component networks in natural rubber. RSC Adv. 2015, 5, 91742–91750. 10.1039/C5RA07428B. [CrossRef] [Google Scholar]
Zhang H.; Zhang L.; Chen X.; Wang Y.; Zhao F.; Luo M.; Liao S. The role of non-rubber components on molecular network of natural rubber during accelerated storage. Polymers 2020, 12, 2880.10.3390/polym12122880. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Kutsukawa R.; Imaizumi R.; Suenaga-Hiromori M.; Takeshita K.; Sakai N.; Misawa S.; Yamamoto M.; Yamaguchi H.; Miyagi-Inoue Y.; Waki T.; Kataoka K.; Nakayama T.; Yamashita S.; Takahashi S. Structure-based engineering of a short-chain cis-prenyltransferase to biosynthesize nonnatural all-cis-polyisoprenoids: molecular mechanisms for primer substrate recognition and ultimate product chain-length determination. FEBS J. 2022, 289, 4602–4621. 10.1111/febs.16392. [Abstract] [CrossRef] [Google Scholar]
Eng A.; Kawahara S.; Tanaka Y. Trans-isofrene units in natural rubber. Rubber Chem. Technol. 1994, 67, 159–168. 10.5254/1.3538662. [CrossRef] [Google Scholar]
Gimenez-Dejoz J.; Tsunoda K.; Fukushima Y.; Numata K. Computational study of the interaction between natural rubber α-terminal groups and l-quebrachitol, one of the major components of natural rubber. Polym. J. (Tokyo, Jpn.) 2022, 54, 229–233. 10.1038/s41428-021-00569-w. [CrossRef] [Google Scholar]
Dixit M.; Taniguchi T.Role of terminal groups of cis-1,4-polyisoprene in the formation of physical cross-linking in NR — All-atom simulation study. In Bulletin of the American Physical Society, APS March Meeting 2023; American Physical Society, 2023; Vol 68
Dixit M.; Taniguchi T. Role of Terminal Groups of cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber. Biomacromolecules 2023, 24, 3589–3602. 10.1021/acs.biomac.3c00355. [Abstract] [CrossRef] [Google Scholar]
Dixit M.; Taniguchi T. Effect of Impurities on the Formation of End-Group Clusters in Natural Rubber: Phenylalanine Dipeptide as an Impurity Protein. Macromolecules 2024, 57, 2588–2608. 10.1021/acs.macromol.3c01833. [CrossRef] [Google Scholar]
Dixit M.; Taniguchi T. Substantial Effect of Terminal Groups in cis-Polyisoprene: A Multiscale Molecular Dynamics Simulation Study. Macromolecules 2022, 55, 9650–9662. 10.1021/acs.macromol.2c01414. [CrossRef] [Google Scholar]
Tangpakdee J.; Tanaka Y. Characterization of sol and gel in hevea natural rubber. Rubber Chem. Technol. 1997, 70, 707–713. 10.5254/1.3538454. [CrossRef] [Google Scholar]
Tanaka Y.; Tarachiwin L. Recent Advances in Structural Characterization of Natural Rubber. Rubber Chem. Technol. 2009, 82, 283–314. 10.5254/1.3548250. [CrossRef] [Google Scholar]
Kitaura T.; Kobayashi M.; Tarachiwin L.; Kum-ourm H.; Matsuura A.; Fushihara K.; Ute K. Characterization of Natural Rubber End Groups Using High-Sensitivity NMR. Macromol. Chem. Phys. 2018, 219, 1700331.10.1002/macp.201700331. [CrossRef] [Google Scholar]
Luo M. C.; Zeng J.; Fu X.; Huang G.; Wu J. Toughening diene elastomers by strong hydrogen bond interactions. Polymer 2016, 106, 21–28. 10.1016/j.polymer.2016.10.056. [CrossRef] [Google Scholar]
Hu B.; Zhou Y.; Luo M. C.; Wei Y. C.; Liu G. X.; Liao S.; Zhao Y. Influence of l -quebrachitol on the properties of centrifuged natural rubber. e-Polym. 2021, 21, 420–427. 10.1515/epoly-2021-0042. [CrossRef] [Google Scholar]
Liu X.-X.; He M.-F.; Luo M.-C.; Wei Y.-C.; Liao S. The role of natural rubber endogenous proteins in promoting the formation of vulcanization networks. e-Polym. 2022, 22, 445–453. 10.1515/epoly-2022-0043. [CrossRef] [Google Scholar]
Wei Y. C.; Liu G. X.; Zhang H. F.; Zhao F.; Luo M. C.; Liao S. Non-rubber components tuning mechanical properties of natural rubber from vulcanization kinetics. Polymer 2019, 183, 121911–121917. 10.1016/j.polymer.2019.121911. [CrossRef] [Google Scholar]
Li S. Q.; Tang M. Z.; Huang C.; Zhang R.; Huang G. S.; Xu Y. X. The Relationship between Pendant Phosphate Groups and Mechanical Properties of Polyisoprene Rubber. Chin. J. Polym. Sci. 2021, 39, 465–473. 10.1007/s10118-021-2497-z. [CrossRef] [Google Scholar]
Zhang R.; Li S.-Q.; Xu R.; Wang C.-C.; Wang Y.; Huang G.; Tang M.; Xu Y.-X. AGGREGATION BEHAVIORS OF PENDANT PHOSPHORYLCHOLINE GROUPS AND THEIR INFLUENCE ON POLYISOPRENE PROPERTIES. Rubber Chem. Technol. 2023, 96, 259–275. 10.5254/rct.23.76948. [CrossRef] [Google Scholar]
Kremer K.; Grest G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. J. Chem. Phys. 1990, 92, 5057–5086. 10.1063/1.458541. [CrossRef] [Google Scholar]
Xu Y.; Hamada Y.; Taniguchi T. Multiscale simulations for polymer melt spinning process using Kremer–Grest CG model and continuous fluid mechanics model. J. Non-Newtonian Fluid Mech. 2024, 325, 105195.10.1016/j.jnnfm.2024.105195. [CrossRef] [Google Scholar]
Moe N. E.; Ediger M. D. Molecular Dynamics Computer Simulation of Polyisoprene Local Dynamics in Dilute Toluene Solution. Macromolecules 1995, 28, 2329–2338. 10.1021/ma00111a028. [CrossRef] [Google Scholar]
Moe N. E.; Ediger M. D. Molecular dynamics computer simulation of local dynamics in polyisoprene melts. Polymer 1996, 37, 1787–1795. 10.1016/0032-3861(96)87294-6. [CrossRef] [Google Scholar]
Moe N. E.; Ediger M. D. Computer simulations of polyisoprene local dynamics in vacuum, solution, and the melt: Are conformational transitions always important?. Macromolecules 1996, 29, 5484–5492. 10.1021/ma960204t. [CrossRef] [Google Scholar]
Faller R.; Reith D. Properties of poly(isoprene): Model building in the melt and in solution. Macromolecules 2003, 36, 5406–5414. 10.1021/ma025877s. [CrossRef] [Google Scholar]
Egami T. Local structure and dynamics of ferroelectric solids. Struct. Bonding (Berlin) 2007, 124, 69–88. 10.1007/430_2006_047. [CrossRef] [Google Scholar]
Diani J.; Gall K. Molecular dynamics simulations of the shape-memory behaviour of polyisoprene. Smart Mater. Struct. 2007, 16, 1575–1583. 10.1088/0964-1726/16/5/011. [CrossRef] [Google Scholar]
Maple J. R.; Dinur U.; Hagler A. T. Derivation of force fields for molecular mechanics and dynamics from ab initio energy surfaces. Proc. Natl. Acad. Sci. U.S.A. 1988, 85, 5350–5354. 10.1073/pnas.85.15.5350. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Deutsch H. P.; Binder K. Interdiffusion and self-diffusion in polymer mixtures: A Monte Carlo study. J. Chem. Phys. 1991, 94, 2294–2304. 10.1063/1.459901. [CrossRef] [Google Scholar]
Ercolessi F.; Adams J. B. Interatomic Potentials from First-Principles Calculations: The Force-Matching Method. Europhys. Lett. 1994, 26, 583–588. 10.1209/0295-5075/26/8/005. [CrossRef] [Google Scholar]
Tschöp W.; Kremer K.; Batoulis J.; Bürger T.; Hahn O. Simulation of polymer melts. I. Coarse-graining procedure for polycarbonates. Acta Polym. 1998, 49, 61–74. 10.1002/(SICI)1521-4044(199802)49:2/3<61::AID-APOL61>3.0.CO;2-V. [CrossRef] [Google Scholar]
Lyubartsev A. P.; Laaksonen A. Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach. Phys. Rev. E 1995, 52, 3730–3737. 10.1103/PhysRevE.52.3730. [Abstract] [CrossRef] [Google Scholar]
Meyer H.; Biermann O.; Faller R.; Reith D.; Müller-Plathe F. Coarse graining of nonbonded inter-particle potentials using automatic simplex optimization to fit structural properties. J. Chem. Phys. 2000, 113, 6264–6275. 10.1063/1.1308542. [CrossRef] [Google Scholar]
Mullinax J. W.; Noid W. G. Reference state for the generalized Yvon-Born-Green theory: Application for coarse-grained model of hydrophobic hydration. J. Chem. Phys. 2010, 133, 124107.10.1063/1.3481574. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Shell M. S. The relative entropy is fundamental to multiscale and inverse thermodynamic problems. J. Chem. Phys. 2008, 129, 144108.10.1063/1.2992060. [Abstract] [CrossRef] [Google Scholar]
Izvekov S.; Voth G. A. A multiscale coarse-graining method for biomolecular systems. J. Phys. Chem. B 2005, 109, 2469–2473. 10.1021/jp044629q. [Abstract] [CrossRef] [Google Scholar]
Uddin M. S.; Ju J. Multiscale modeling of a natural rubber: Bridging a coarse-grained molecular model to the rubber network theory. Polymer 2016, 101, 34–47. 10.1016/j.polymer.2016.08.037. [CrossRef] [Google Scholar]
Li Y.; Kröger M.; Liu W. K. Primitive chain network study on uncrosslinked and crosslinked cis-polyisoprene polymers. Polymer 2011, 52, 5867–5878. 10.1016/j.polymer.2011.10.044. [CrossRef] [Google Scholar]
Kawaguchi N.; Okubo K.; Aida K. Molecular cloning and expression of corticotropin-releasing hormone and urotensin I in the medaka, Oryzias latipes. Fish. Sci. 2002, 68, 1281–1282. 10.2331/fishsci.68.sup2_1281. [CrossRef] [Google Scholar]
Pandey Y. N.; Brayton A.; Burkhart C.; Papakonstantopoulos G. J.; Doxastakis M. Multiscale modeling of polyisoprene on graphite. J. Chem. Phys. 2014, 140, 054908.10.1063/1.4863918. [Abstract] [CrossRef] [Google Scholar]
Moe N. E.; Ediger M. D. Calculation of the coherent dynamic structure factor of polyisoprene from molecular dynamics simulations. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 59, 623–630. 10.1103/PhysRevE.59.623. [CrossRef] [Google Scholar]
Li Y.; Mattice W. L. Atom-Based Modeling of Amorphous 1,4-cis-Polybutadiene. Macromolecules 1992, 25, 4942–4947. 10.1021/ma00045a020. [CrossRef] [Google Scholar]
Ohkuma T.; Kremer K. A composition transferable and time-scale consistent coarse-grained model for cis-polyisoprene and vinyl-polybutadiene oligomeric blends. JPhys Mater. 2020, 3, 034007.10.1088/2515-7639/ab906b. [CrossRef] [Google Scholar]
Ohkuma T.; Kremer K. Comparison of two coarse-grained models of cis-polyisoprene with and without pressure correction. Polymer 2017, 130, 88–101. 10.1016/j.polymer.2017.09.062. [CrossRef] [Google Scholar]
Shukla R.; Ray D.; Sarkar K.; Kumar Dixit M.; Prasad Bhattacharyya S. Flying onto global minima on potential energy surfaces: A swarm intelligence guided route to molecular electronic structure. Int. J. Quantum Chem. 2016, 117, e2532810.1002/qua.25328. [CrossRef] [Google Scholar]
Ghosh S.; Dixit M.; Bhattacharyya S.; Tembe B. Franck-condon factors for diatomics: Insights and analysis using the Fourier Grid Hamiltonian method. J. Chem. Educ. 2013, 90, 1463–1471. 10.1021/ed4002199. [CrossRef] [Google Scholar]
Siddique A.; Dixit M.; Tembe B. Molecular dynamics simulations of Ca²⁺ Cl^– ion pair in polar mixtures of acetone and water: Solvation and dynamical studies. Chem. Phys. Lett. 2016, 662, 306–316. 10.1016/j.cplett.2016.09.061. [CrossRef] [Google Scholar]
Choi Y. K.; Park S. J.; Park S.; Kim S.; Kern N. R.; Lee J.; Im W. CHARMM-GUI Polymer Builder for Modeling and Simulation of Synthetic Polymers. J. Chem. Theory Comput. 2021, 17, 2431–2443. 10.1021/acs.jctc.1c00169. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Lehman N.; Tuljittraporn A.; Songtipya L.; Uthaipan N.; Sengloyluan K.; Johns J.; Nakaramontri Y.; Kalkornsurapranee E. Influence of Non-Rubber Components on the Properties of Unvulcanized Natural Rubber from Different Clones. Polymers 2022, 14, 1759.10.3390/polym14091759. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Nun-anan P.; Wisunthorn S.; Pichaiyut S.; Nathaworn C. D.; Nakason C. Influence of nonrubber components on properties of unvulcanized natural rubber. Polym. Adv. Technol. 2020, 31, 44–59. 10.1002/pat.4746. [CrossRef] [Google Scholar]
Nun-anan P.; Wisunthorn S.; Pichaiyut S.; Vennemann N.; Nakason C. Novel approach to determine non-rubber content in Hevea brasiliensis: Influence of clone variation on properties of un-vulcanized natural rubber. Ind. Crops Prod. 2018, 118, 38–47. 10.1016/j.indcrop.2018.03.011. [CrossRef] [Google Scholar]
Mooibroek H.; Cornish K. Alternative sources of natural rubber. Appl. Microbiol. Biotechnol. 2000, 53, 355–365. 10.1007/s002530051627. [Abstract] [CrossRef] [Google Scholar]
Lockwood G. GAS CHROMATOGRAPHY | Terpenoids. Reference Module in Chemistry, Molecular Sciences and Chemical Engineering 2013, 10.1016/b978-0-12-409547-2.04768-5. [CrossRef] [Google Scholar]
Abraham M. J.; Murtola T.; Schulz R.; Páll S.; Smith J. C.; Hess B.; Lindahl E. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1–2, 19–25. 10.1016/j.softx.2015.06.001. [CrossRef] [Google Scholar]
Allouche A.-r. Gabedit—A graphical user interface for computational chemistry softwares. J. Comput. Chem. 2011, 32, 174–182. 10.1002/jcc.21600. [Abstract] [CrossRef] [Google Scholar]
Vanommeslaeghe K.; Hatcher E.; Acharya C.; Kundu S.; Zhong S.; Shim J.; Darian E.; Guvench O.; Lopes P.; Vorobyov I.; Mackerell A. D. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671–690. 10.1002/jcc.21367. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Parrinello M.; Rahman A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182–7190. 10.1063/1.328693. [CrossRef] [Google Scholar]
Hess B.; Bekker H.; Berendsen H. J.; Fraaije J. G. LINCS: A Linear Constraint Solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463–1472. 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H. [CrossRef] [Google Scholar]
Darden T.; York D.; Pedersen L. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092. 10.1063/1.464397. [CrossRef] [Google Scholar]
Fetters L.; Lohse D.; Colby R.Physical Properties of Polymers Handbook; Springer, 2007; pp 447–454. [Google Scholar]
Williams G.; Watts D. C. Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function. Trans. Faraday Soc. 1970, 66, 80–85. 10.1039/tf9706600080. [CrossRef] [Google Scholar]
Doi M.; Edwards S.; Edwards S.The Theory of Polymer Dynamics; International Series of Monographs on Physics; Clarendon Press, 1988; . [Google Scholar]
Verdier P. H. Monte carlo studies of lattice-model polymer chains. i. correlation functions in the statistical-bead model. J. Chem. Phys. 1966, 45, 2118–2121. 10.1063/1.1727896. [CrossRef] [Google Scholar]
Kopf A.; Dünweg B.; Paul W. Dynamics of polymer ”isotope” mixtures: Molecular dynamics simulation and Rouse model analysis. J. Chem. Phys. 1997, 107, 6945–6955. 10.1063/1.474934. [CrossRef] [Google Scholar]
Harmandaris V. A.; Mavrantzas V. G.; Theodorou D. N. Atomistic molecular dynamics simulation of polydisperse linear polyethylene melts. Macromolecules 1998, 31, 7934–7943. 10.1021/ma980698p. [CrossRef] [Google Scholar]
Tsolou G.; Harmandaris V. A.; Mavrantzas V. G. Molecular dynamics simulation of temperature and pressure effects on the intermediate length scale dynamics and zero shear rate viscosity of cis-1,4-polybutadiene: Rouse mode analysis and dynamic structure factor spectra. J. Non-Newtonian Fluid Mech. 2008, 152, 184–194. 10.1016/j.jnnfm.2007.10.011. [CrossRef] [Google Scholar]
Vandoolaeghe W. L.; Terentjev E. M. Constrained Rouse model of rubber viscoelasticity. J. Chem. Phys. 2005, 123, 034902.10.1063/1.1955445. [Abstract] [CrossRef] [Google Scholar]
Janosi L.; Li Z.; Hancock J. F.; Gorfe A. A. Organization, dynamics, and segregation of Ras nanoclusters in membrane domains. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 8097–8102. 10.1073/pnas.1200773109. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Li Z.; Dormidontova E. E. Kinetics of diblock copolymer micellization by dissipative particle dynamics. Macromolecules 2010, 43, 3521–3531. 10.1021/ma902860j. [CrossRef] [Google Scholar]
Marrink S. J.; Tieleman D. P.; Mark A. E. Molecular dynamics simulation of the kinetics of spontaneous micelle formation. J. Phys. Chem. B 2000, 104, 12165–12173. 10.1021/jp001898h. [CrossRef] [Google Scholar]
Dixit M.; Chatterjee A.; Tembe B. Salting-out of methane in the aqueous solutions of urea and sarcosine. J. Chem. Sci. 2016, 128, 599–612. 10.1007/s12039-016-1052-x. [CrossRef] [Google Scholar]
Dixit M.; Hajari T.; Tembe B. The effect of urea and taurine osmolytes on hydrophobic association and solvation of methane and neopentane molecules. J. Mol. Liq. 2016, 223, 660–671. 10.1016/j.molliq.2016.08.079. [CrossRef] [Google Scholar]
Dixit M.; Hajari T.; Tembe B. Solvation structures of sodium halides in dimethyl sulfoxide (DMSO) – – methanol (MeOH) mixtures. Mol. Simul. 2017, 43, 154–168. 10.1080/08927022.2016.1241396. [CrossRef] [Google Scholar]
Dixit M.; Hajari T.; Meti M. D.; Srivastava S.; Srivastava A.; Daniel J. Ionic Pairing and Selective Solvation of Butylmethylimidazolium Chloride Ion Pairs in DMSO–Water Mixtures: A Comprehensive Examination via Molecular Dynamics Simulations and Potentials of Mean Force Analysis. J. Phys. Chem. B 2024, 128, 2168–2180. 10.1021/acs.jpcb.3c06876. [Abstract] [CrossRef] [Google Scholar]
Dixit M.; Siddique A.; Tembe B. Salting-Out of Methane in the Aqueous Solutions of Urea and Glycine-Betaine. J. Phys. Chem. B 2015, 119, 10941–10953. 10.1021/acs.jpcb.5b00556. [Abstract] [CrossRef] [Google Scholar]
Sarkar A.; Dixit M.; Tembe B. Solvation structures of lithium halides in methanol-water mixtures. Chem. Phys. 2015, 447, 76–85. 10.1016/j.chemphys.2014.11.019. [CrossRef] [Google Scholar]
Meti M.; Dixit M.; Hajari T.; Tembe B. Ion pairing and preferential solvation of butylmethylimidazolium chloride ion pair in water-ethanol mixtures by using molecular dynamics simulations. Chem. Phys. Lett. 2019, 720, 107–112. 10.1016/j.cplett.2019.02.010. [CrossRef] [Google Scholar]
Hajari T.; Dixit M.; Yadav H. Hydrophobic association and solvation of neopentane in urea, TMAO and urea-TMAO solutions. Phys. Chem. Chem. Phys. 2022, 24, 6941–6957. 10.1039/D1CP05321C. [Abstract] [CrossRef] [Google Scholar]
Ghosh S.; Dixit M.; Chakrabarti R. Thermodynamics of site-specific small molecular ion interactions with DNA duplex: A molecular dynamics study. Mol. Simul. 2016, 42, 715–724. 10.1080/08927022.2015.1085123. [CrossRef] [Google Scholar]
Meti M.; Dixit M.; Tembe B. Salting-in of neopentane in the aqueous solutions of urea and glycine-betaine. Mol. Simul. 2018, 44, 677–687. 10.1080/08927022.2018.1431834. [CrossRef] [Google Scholar]
Kumar A.; Mahato J.; Dixit M.; Patwari G. Progressive Hydrophobicity of Fluorobenzenes. J. Phys. Chem. B 2019, 123, 10083–10088. 10.1021/acs.jpcb.9b08057. [Abstract] [CrossRef] [Google Scholar]
Chatterjee A.; Dixit M.; Tembe B. Solvation Structures and Dynamics of the Magnesium Chloride (Mg²⁺ - Cl^–) Ion Pair in Water-Ethanol Mixtures. J. Phys. Chem. A 2013, 117, 8703–8709. 10.1021/jp4031706. [Abstract] [CrossRef] [Google Scholar]
Dixit M.; Tembe B. Potentials of mean force of sodium chloride ion pair in dimethyl sulfoxide-methanol mixtures. J. Mol. Liq. 2013, 178, 78–83. 10.1016/j.molliq.2012.09.026. [CrossRef] [Google Scholar]
Siddique A.; Dixit M.; Tembe B. Solvation structure and dynamics of potassium chloride ion pair in dimethyl sulfoxide-water mixtures. J. Mol. Liq. 2013, 188, 5–12. 10.1016/j.molliq.2013.09.004. [CrossRef] [Google Scholar]
Jain A.; Dixit M.; Tembe B. Preferential solvation and association constants of 2- Exo and 2- Endo norbornyl chlorides in water-acetone mixtures. Mol. Simul. 2014, 40, 987–995. 10.1080/08927022.2013.832246. [CrossRef] [Google Scholar]
Dixit M.; Lazaridis T. Free energy of hydrophilic and hydrophobic pores in lipid bilayers by free energy perturbation of a restraint. J. Chem. Phys. 2020, 153, 054101.10.1063/5.0016682. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Justice M. C.; Justice J. C. Ionic interactions in solutions. I. The association concepts and the McMillan-Mayer theory. J. Solution Chem. 1976, 5, 543–561. 10.1007/BF00647377. [CrossRef] [Google Scholar]
Liu P.; Harder E.; Berne B. J. On the calculation of diffusion coefficients in confined fluids and interfaces with an application to the liquid-vapor interface of water. J. Phys. Chem. B 2004, 108, 6595–6602. 10.1021/jp0375057. [CrossRef] [Google Scholar]
Araya-Secchi R.; Perez-Acle T.; Kang S. G.; Huynh T.; Bernardin A.; Escalona Y.; Garate J. A.; Martínez A.; García I.; Sáez J.; Zhou R. Characterization of a novel water pocket inside the human Cx26 hemichannel structure. Biophys. J. 2014, 107, 599–612. 10.1016/j.bpj.2014.05.037. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
Gowers R.; Linke M.; Barnoud J.; Reddy T.; Melo M.; Seyler S.; Domański J.; Dotson D.; Buchoux S.; Kenney I.; Beckstein O.MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. In Proceedings of the 15th Python in Science Conference, 2016; pp 98–105.
Paul W.; Smith G. D.; Yoon D. Y. Static and dynamic properties of a n-C100H202 melt from molecular dynamics simulations. Macromolecules 1997, 30, 7772–7780. 10.1021/ma971184d. [CrossRef] [Google Scholar]
Paul W.; Smith G. D.; Yoon D. Y.; Farago B.; Rathgeber S.; Zirkel A.; Willner L.; Richter D. Chain motion in an unentangled polyethylene melt: A critical test of the rouse model by molecular dynamics simulations and neutron spin echo spectroscopy. Phys. Rev. Lett. 1998, 80, 2346–2349. 10.1103/PhysRevLett.80.2346. [CrossRef] [Google Scholar]
Paul W.; Binder K.; Heermann D. W.; Kremer K. Dynamics of polymer solutions and melts. Reptation predictions and scaling of relaxation times. J. Chem. Phys. 1991, 95, 7726–7740. 10.1063/1.461346. [CrossRef] [Google Scholar]
Smith G. D.; Paul W. United Atom Force Field for Molecular Dynamics Simulations of 1,4-Polybutadiene Based on Quantum Chemistry Calculations on Model Molecules. J. Phys. Chem. A 1998, 102, 1200–1208. 10.1021/jp9730858. [CrossRef] [Google Scholar]
Smith G. D.; Paul W.; Monkenbusch M.; Richter D. On the non-Gaussianity of chain motion in unentangled polymer melts. J. Chem. Phys. 2001, 114, 4285–4288. 10.1063/1.1348032. [CrossRef] [Google Scholar]
Rapaport D. C. Hydrogen bonds in water network organization and lifetimes. Mol. Phys. 1983, 50, 1151–1162. 10.1080/00268978300102931. [CrossRef] [Google Scholar]
Skarmoutsos I.; Guardia E.; Samios J. Hydrogen bond, electron donor-acceptor dimer, and residence dynamics in supercritical CO2-ethanol mixtures and the effect of hydrogen bonding on single reorientational and translational dynamics: A molecular dynamics simulation study. J. Chem. Phys. 2010, 133, 014504.10.1063/1.3449142. [Abstract] [CrossRef] [Google Scholar]
Skarmoutsos I.; Guardia E. Local Structural Effects and Related Dynamics in Supercritical Ethanol. 2. Hydrogen-Bonding Network and Its Effect on Single Reorientational Dynamics. J. Phys. Chem. B 2009, 113, 8898–8910. 10.1021/jp901489c. [Abstract] [CrossRef] [Google Scholar]
Stillinger F. H.. Advances in Chemical Physics; Wiley, 1975; Vol. 31. [Google Scholar]
Gowers R. J.; Carbone P. A multiscale approach to model hydrogen bonding: The case of polyamide. J. Chem. Phys. 2015, 142, 224907.10.1063/1.4922445. [Abstract] [CrossRef] [Google Scholar]

Articles from ACS Polymers Au are provided here courtesy of American Chemical Society

Full text links

Read article at publisher's site: https://doi.org/10.1021/acspolymersau.4c00019

Search life-sciences literature (45,103,589 articles, preprints and more)

Exploring the Role of Hydroxy- and Phosphate-Terminated cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber: Insights from Molecular Dynamics Simulations.

Affiliations

Authors

ORCIDs linked to this article

Abstract

Free full text

Exploring the Role of Hydroxy- and Phosphate-Terminated cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber: Insights from Molecular Dynamics Simulations

Associated Data

Abstract

1. Introduction

Table 1

2. Methodology and Computational Discriptions

2.1. All-Atom Molecular Dynamics Protocol

3. Results

3.1. Equilibration of Polyisoprene Chains

Table 2

3.2. Rouse Mode Analysis

Table 3

3.3. Structure Aggregates of Hydroxylated- and Phosphorylated cis-1,4-Polyisoprene Chains

3.4. Survival Time of Terminal Groups in the Coordination Shell

Table 4

3.5. Self-Diffusion Behavior

3.6. Characterization of Polar Physical Junction Points: Insights into Structural Terminal Groups

Table 5

4. Discussion and Concluding Remarks

Acknowledgments

Supporting Information Available

Notes

Notes

Supplementary Material

References

Full text links

Similar Articles

Role of Terminal Groups of cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber.

Structural Analysis of the Terminal Groups in Commercial Hevea Natural Rubber by 2D-NMR with DOSY Filters and Multiple-WET Methods Using Ultrahigh-Field NMR.

Nanocomposite of Fullerenes and Natural Rubbers: MARTINI Force Field Molecular Dynamics Simulations.

Structural characterization of natural rubber based on recent evidence from selective enzymatic treatments.

Partnerships & funding

Similar Articles

Role of Terminal Groups of cis-1,4-Polyisoprene Chains in the Formation of Physical Junction Points in Natural Rubber.
Biomacromolecules, 24(8):3589-3602, 01 Aug 2023

Structural Analysis of the Terminal Groups in Commercial Hevea Natural Rubber by 2D-NMR with DOSY Filters and Multiple-WET Methods Using Ultrahigh-Field NMR.
Biomacromolecules, 20(3):1394-1400, 26 Feb 2019

Nanocomposite of Fullerenes and Natural Rubbers: MARTINI Force Field Molecular Dynamics Simulations.

Structural characterization of natural rubber based on recent evidence from selective enzymatic treatments.
J Biosci Bioeng, 103(4):287-292, 01 Apr 2007