Plasma and Fusion Research: Regular Articles

Volume 14, 3401106 (2019)

Structural Changes in Tritium-Substituted Polymeric Materials by
Beta Decays: A Molecular Dynamics Study∗)
Haolun LI, Susumu FUJIWARA, Hiroaki NAKAMURA1,2) , Tomoko MIZUGUCHI,Takuo YASUNAGA3) , Takao OTSUKA4) , Takahiro KENMOTSU5) , Yuji HATANO6)and Shinji SAITO7)
Graduate School of Science and Technology, Kyoto Institute of Technology, Kyoto 606-8585, Japan1)
Graduate School of Engineering, Nagoya University, Nagoya 464-8601, Japan
2)
National Institute for Fusion Science, National Institutes of Natural Sciences, Toki 509-5292, Japan3)
Kyushu Institute of Technology, Iizuka 820-8502, Japan
4)
RIKEN, Wako 351-0198, Japan
5)
Doshisha University, Kyotanabe 610-0394, Japan
6)
University of Toyama, Toyama 930-8555, Japan
7)
Institute for Molecular Science, National Institutes of Natural Sciences, Okazaki 444-8585, Japan(Received 10 January 2019 / Accepted 10 April 2019)
The molecular mechanism through which how beta decays in tritium-substituted species damage DNA andpolymeric materials is still unknown. Molecular dynamics simulations of hydrogen-removed polyethylene wereperformed to predict the structural change of the polyethylene chain after the substituted tritium decays. We
calculated the potential energy, the global orientational order parameter, and the average number of consecutive
trans bonds. The results are that, the greater the number of removed hydrogen atoms, the higher the potential
energy and the lower the value of the global orientational order parameter and the average number of consecutive
trans bonds. Thus, after losing hydrogen, polyethylene becomes poorer in terms of both thermal and structuralstabilities.
c 2019 The Japan Society of Plasma Science and Nuclear Fusion ResearchKeywords: polyethylene, molecular dynamics simulation, beta decay, structural change, tritium
DOI: 10.1585/pfr.14.3401106

1. Introduction

by beta decays of substituted tritium.

Following the Fukushima nuclear accident in 2011,a large amount of tritiated water was discharged, whichcaused high levels of pollution to the surrounding land andthe sea. Tritiated water is a type of water which contains tritium (T or 3 H), a radioactive isotope of hydrogen. With a half-life of 12.323 years, tritium decays tohelium-3 (3 He) by beta decay with emissions of a betaelectron (beta-ray) and an antineutrino [1]. In particular,
on macromolecules, such as DNA and polymeric materials, beta-rays may cause damage, and if substituted tritium decays to helium-3, chemical bonds break, destabilizing the structure. However, the molecular mechanismof this damage caused by tritium is still an unsolved problem. Therefore, our aim is to predict the structural changeof tritium-substituted macromolecules, such as DNA andpolymeric materials, by beta decays of tritium to helium-3
using molecular dynamics (MD) simulations.
In this paper, MD simulations of polyethylene were
performed to study the structural change of polyethylene
author’s e-mail: m8672034@edu.kit.ac.jp∗) This article is based on the presentation at the 27th International Toki
Conference (ITC27) & the 13th Asia Pacific Plasma Theory Conference(APPTC2018).

2. Model and Simulation
The polyethylene model is made up of a linear chainof 2998 CH2 groups with two terminal CH3 groups alltreated as united atoms, with masses of 14 and 15 g/mol,respectively. After decay of substituted tritium to helium3, the C–T bonds break, as the inert helium-3 is freed. To
imitate the decay eﬀect on the polymer chain, hydrogenatoms are removed randomly from the polyethylene. The
remaining -CH- united atoms, with a mass of 13 g/mol, areconsidered overlapping sp2 hybrid orbitals [2]. The ratioof the number of removed hydrogen atoms to the numberof carbon atoms, fH , is set at fH = 0, 0.001, 0.01 and 0.1.
The DREIDING potential [3] is the generic force fieldused in this research, which is useful for predicting structures and dynamics of polymers. As atom types of theDREIDING force field, a label such as “C_32” is used tostand for each type of united atom. The first two characters correspond to the chemical elements (e.g., “He” forhelium), where an underscore is added for elements with
only one letter (e.g., “C_” for carbon). The third charactersuggests the hybridization: 1 for linear (sp1 ), 2 for trigonal

3401106-1

c 2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 14, 3401106 (2019)
Table 3 Parameters used for the torsion interaction.

Table 1 Equilibrium bond distance R0I J , where “C_3” includes both “C_33” and “C_32”, and “C_2” represents“C_21”.

Table 4 Parameters used for the Lennard–Jones potential.
Table 2 Equilibrium angle θ0J , where “C_3” includes both
“C_33” and “C_32”, and “C_2” represents “C_21”.

(sp2 ), and 3 for tetrahedral (sp3 ). The fourth character indicates the number of implicit hydrogen atoms. As part ofa united atom, hydrogen atoms are not to be calculated independently. Therefore, there are three atom types used inthis paper: “C_33” for the CH3 group, “C_32” for the CH2
group, and “C_21” for the -CH- group, which is generated
by removing a hydrogen atom from the CH2 group.
The potential energy (1) includes bonded interactions(Eval ) and nonbonded interactions (Enb ):Epot = Eval + Enb .

(1)

In this paper, Eval (2) consists of the sum of bond stretch(EB ), bond-angle bend (EA ), and dihedral angle torsion(ET ),
Eval = EB + EA + ET ,(2)

whereas Enb is the Lennard–Jones potential:
Enb = ELJ .

(3)

The bond stretch interaction (EB ) is given by the expression (4)EB =1
KI J (RI J − R0I J )2 ,2(4)where RI J is the bond distance between atom I and atom J.
The force constant is KI J = 700 (kcal/mol)/Å2 . R0I J is theequilibrium bond distance, and Table 1 shows the values ofR0I J .
The angle bend interaction is calculated using Eq. (5):EA =1KI JK (θI JK − θ0J )2 ,2(5)where θI JK is the angle between bonds I J and JK, and theforce constant is KI JK = 100 (kcal/mol)/rad2 for the angle
bend. θ0J represents the equilibrium angle. The values of
θ0J are given in Table 2.
The torsion interaction (6) takes the form

ET =1
V JK {1 − cos[n JK (ϕ − ϕ0JK )]},2(6)

where ϕ is the dihedral angle between the planes I JK andJKL, V JK is taken as the energy barrier to rotation, n JK isthe periodicity, and ϕ0JK is the equilibrium angle. As V JK ,n JK , and ϕ0JK are independent of atom I and atom L, thevalues are as shown in Table 3.
For the Lennard–Jones potential as a nonbonded interaction, the common Eq. (7) is used: 12  6 σσ−,(7)
ELJ (r) = 4εrr
with a cutoﬀ of 10.5 Å, where r is the distance between
two united atoms, ε is the depth of the potential well andσ is the finite distance at which the interparticle potential
is zero. Table 4 shows the values of ε and σ.
By using the velocity version of the Verlet algorithm[4], the equations of motion for all united atoms are solved
numerically. To keep the temperature of the system constant, the Nosé–Hoover method is applied [5–7]. The integration time step is 1 fs.
Three diﬀerent simulation runs are performed for eachfH as follows. First, for preparation, the isolated polyethylene is cooled from 800 to 300 K, at a rate of 50 K/ns
(10,000,000 time steps or 10 ns in total) in a vacuum
simulation box, forming a folded, orientationally orderedstructure, in agreement with S. Fujiwara and T. Sato [8].
Secondly, hydrogen atoms are randomly removed fromthe polyethylene model for three simulation runs. For
fH = 0, the simulations are run for an extra 1,000,000 and
2,000,000 time steps (1 ns and 2 ns, respectively) at 300 K
as a preparation for the second and third simulation runs,respectively. After that, the model is heated stepwise to550 K, at a rate of 10 K/ns (1,000,000 time steps at each3401106-2

Plasma and Fusion Research: Regular Articles

Volume 14, 3401106 (2019)

temperature, 25,000,000 time steps or 25 ns in total), andthe model still maintains a folded, orientationally orderedstructure at 550 K. Then, the temperature of the simulation
box undergoes a second stepwise increase, to 800 K, at thelow rate of 2 K/ns (1,000,000 time steps at each temperature, 125,000,000 time steps or 125 ns in total). All thesimulations were run using LAMMPS [9, 10].

3. Results and Discussion
3.1

Structural changes

Structures were visualized using the VMD software[11, 12]. Structural changes in all models are observed,
from an orientationally ordered structure to a disorderedone. Figures 1 and 2 show the structure of the models forone of the simulation runs before and after the structuralchange at a certain temperature of 650 and 700 K, respectively. In the latter, the main chains are shown in blue lines,
whereas -CH- groups are shown in red balls.

3.2

Potential energy

As mentioned in Section 2, the potential energy is calculated using Eqs. (1) - (7), and Fig. 3 shows the average
curves obtained for three diﬀerent simulation runs.
It can be seen from Fig. 3 that curves for fH = 0,
fH = 0.001 and fH = 0.01 are similar, although Epot atfH = 0.1 is the largest among the range of fH tested. Epotfor fH = 0.1 increases sharply from about 650 to 690 K,
whereas the sharp increase of Epot for other fH ’s starts atabout 650 K and remains up to about 720 K. Above 720 K,
all Epot curves converge to a linearity, in a similar fashion. This observation can be interpreted as follows. We
find that during heating, the values of EB for diﬀerent values of fH are almost the same, as well as the values ofEA and ET . However, the values of ELJ for diﬀerent values of fH are diﬀerent. In the temperature range wherethe structure changes, the diﬀerence of ELJ becomes maximum, whereas over 720 K, ELJ becomes almost the same,once the obtained structure is disordered. As a result, theEpot (= EB + EA + ET + ELJ ) over 720 K becomes identical.

3.3 Global orientational order parameter
The global orientational order parameter (P2 ) is calculated using expression (8) to analyze the degree of orientational order in the polyethylene chain as3 cos2 Ψ − 1,
(8)
P2 =2
bond
where Ψ is the angle between two chord vectors and· · · bond denotes the average among all pairs of chord vectors. Here the chord vector is defined as the vector fromthe center of one bond to the center of another adjacentbond along the polymer chain. When all chord vectors areparallel, P2 = 1, and when all chord vectors are randomly
oriented, P2 = 0. Figure 4 shows the dependence of P2 ontemperature T for all fH ’s.
For fH = 0, P2 decreases most quickly at the increasing temperatures from about 665 to 695 K, whereas
for fH  0, P2 drops from 665, 660, 640 K to about 690,690, 670 K for fH = 0.001, fH = 0.01, and fH = 0.1, re-Fig. 2 Side views of polyethylene at 700 K for (a) fH = 0, (b)fH = 0.001, (c) fH = 0.01, (d) fH = 0.1.

Fig. 1 Side views of polyethylene at 650 K for (a) fH = 0, (b)fH = 0.001, (c) fH = 0.01, (d) fH = 0.1.

Fig. 3 Temperature dependence of the potential energy Epot for(a) fH = 0, (b) fH = 0.001, (c) fH = 0.01, (d) fH = 0.1.

3401106-3

Plasma and Fusion Research: Regular Articles

Volume 14, 3401106 (2019)

Fig. 5 The average size of the trans segments ntr vs inverse temperature 1/T for (a) fH = 0, (b) fH = 0.001, (c) fH = 0.01,(d) fH = 0.1. The dashed line is the line fitted for hightemperature region by Eq. (9) with Δε = 2.0 kcal/mol.

Fig. 4 Temperature dependence of P2 for (a) fH = 0, (b) fH =0.001, (c) fH = 0.01, (d) fH = 0.1.

fH = 0.1. The curves are well fitted by Eq. (9) at the hightemperature range over 700 K, with Δε = 2.0 kcal/mol,which is the value of V JK for torsion interaction between
C_32 and C_32 in Eq. (6).

4. Conclusions
spectively. It can be found from the temperature range inwhich such curves drop that the greater the value of fH is,the ordered structure breaks down at lower temperatures.
For all values of fH , there is a clear P2 residual after themain drop, which does not hold for long, disappearing at
about 720 K.

3.4

Conformational change

Moreover, all local structures, consisting of several
united atoms, were checked to investigate conformationalchanges. As the values of the dihedral angle ϕ were examined, the trans state and the gauche state were defined by|ϕ| ≤ π/3 and |ϕ| > π/3, respectively [8].
The average size of trans segments, i.e., the average
number of consecutive trans bonds ntr , was calculated andare shown in Fig. 5. The curves of ntr are the average ofthree simulation runs. ntr corresponds to the persistence
length, which fits to the following relation (9) at high temperatures:Δεntr ∝ exp,(9)kB T
where Δε is the torsional energy barrier between the transand the gauche conformations and kB is the Boltzmann
constant.
As Fig. 5 shows, the values of ntr for fH = 0, fH =0.001, and fH = 0.01 are at the same level, whereas thevalues of ntr for fH = 0.1 are smaller. The curves increase
sharply with decreasing temperature at about 720 K for
fH = 0, fH = 0.001, and fH = 0.01 and at about 690 K for

In this research, the structural changes of tritiumsubstituted polyethylene with three types of united atoms
by beta decays of tritium to helium-3 have been studied
during the gradual stepwise heating by MD simulations.
All models kept a folded, orientationally orderedstructure at first, which became gradually disordered uponheating. As fH becomes larger, it is found that:(i) The visualized structure is more disordered at thesame temperature what is also reflected in the lesser
value of P2 for a greater fH .
(ii) The change of ordered structure starts and ends at alower temperature for a greater value of fH .
(iii) In Figs. 3 - 5, the curves change at the same temperature range for each fH , respectively. Such range canbe regarded as the range where the structural changes
are occurring for each fH .
The results suggest that both the thermal and structural
stabilities of polyethylene can be decreased by removinghydrogen atoms from the chain. The greater the numberof removed hydrogen atoms, the clearer the eﬀect. We alsofind that the curves of physical quantities (Epot , P2 , and ntr )of each simulation run are slightly diﬀerent in shape fromother simulation runs, suggesting that the positions of theremoved hydrogen atoms have an eﬀect on the structure.
In a future work, ReaxFF (Reactive Force Field) MD
simulation studies will be performed to treat all hydrogenatoms and residual electrons properly.

3401106-4

Plasma and Fusion Research: Regular Articles

Volume 14, 3401106 (2019)

Acknowledgments
This study was partially supported by the Joint Studies Program (2017 - 2019) of the Institute for Molecular
Science, the Joint Research by the National Institutes ofNatural Sciences (NINS), JSPS KAKENHI Grant Number
JP15K05244, and the NIFS Collaborative Research Program (NIFS17KNTS050).
[1] T. Tanabe (ed.), Tritium: Fuel of Fusion Reactors (SpringerJapan, Japan, 2017).
[2] S. Fujiwara, H. Nakamura, H. Li, H. Miyanishi, T.
Mizuguchi, T. Yasunaga, T. Otsuka, Y. Hatano and S. Saito,Adv. Simulat. Sci. Eng. 6, 94 (2019).
[3] S.L. Mayo, B.D. Olafson and W.A. Goddard III, J. Phys.
Chem. 94, 8897 (1990).
[4] W.C. Swope, H.C. Andersen, P.H. Berens and K.R. Wilson,J. Chem. Phys. 76, 637 (1982).
[5] S. Nosé, Mol. Phys. 52, 255 (1984).
[6] S. Nosé, J. Chem. Phys. 81, 511 (1984).
[7] W.G. Hoover, Phys. Rev. A 31, 1695 (1985).
[8] S. Fujiwara and T. Sato, J. Chem. Phys. 114, 6455 (2001).
[9] https://lammps.sandia.gov/[10] S. Plimpton, J. Comp. Phys. 117, 1 (1995).
[11] https://www.ks.uiuc.edu/Research/vmd/[12] W. Humphrey, A. Dalke and K. Schulten, J. Molec. Graphics 14, 33 (1996).

3401106-5