Molecular Simulation Studies on the Interactions of 2,4,6

Subscriber access provided by BUFFALO STATE

B: Biomaterials and Membranes

Molecular Simulation Studies on the Interactions of 2,4,6Trinitrotoluene and Its Metabolites with Lipid Membranes Hong Yang, Huarong Li, Liu Liu, Yang Zhou, and Xinping Long J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b03033 • Publication Date (Web): 08 Jul 2019 Downloaded from pubs.acs.org on July 18, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Simulation Studies on the Interactions of 2,4,6-Trinitrotoluene and Its Metabolites with Lipid Membranes Hong Yang,†,‡,‖ Huarong Li,‡,‖ Yang Zhou,‡,* Xinping Long‡

† School

of Material Science and Engineering, Tsinghua University, Beijing 100084, China;

‡ Institute

of Chemical Materials, China Academy of Engineering and Physics, Mianyang 621900, China.

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT The penetration of chemicals into biological membranes is a key factor in the determination of their possible effects on organisms, this complex process mainly involves in the interaction between chemicals and membranes. Here, we reported the interaction between a highly toxic class of explosives (TNT and its metabolites) and lipid membranes using molecular dynamics simulations. We calculated the permeability coefficient, transmembrane time and liposome-water partition coefficient by integrating free energy curves for all species. The results showed that TNT had the lower transmembrane capacity than its metabolites. Based on the liposome-water partition coefficient, we demonstrated that the membrane affinity of TNT is larger than its diamino metabolites, but less than its monoamino metabolites. This result can qualitatively explain the difference of bioconcentration factors in experiments. The accumulation of TNT and metabolites in membranes can change the membrane structure, such as the area per lipid, the thickness of lipid bilayers and the order of lipid tails, further the penetration of water. All these are closely related to the interactions (mainly hydrogen bonds) of TNT and metabolites with lipid and water molecules. This work has a certain significance for understanding the toxicity of TNT and its metabolites.

2


Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


 INTRODUCTION The compound 2, 4, 6-trinitrotoluene (TNT) is a widely used explosive throughout the world, and also a highly toxic xenobiotic. TNT has been rated as a class C carcinogen by the Environment Protection Agency, and can cause the methemoglobinemia, dermatitis, hepatitis, anemia, cataracts, and jaundice for humans.1 Importantly, some diseases have not been alleviated over time, such as cataracts (see the epidemiological survey from China in Table S1). Therefore, the toxicity mechanism of TNT and its metabolites have been always paid more attention.2-5 In the different organisms, TNT can be transformed to different metabolites, including 2-amino-4, 6-dinitrotoluene (2A), 4-amino-2, 6-dinitrotoluene (4A), 2, 4-diamino-6-nitrotoluene (24DA), 2, 6-diamino-4-nitrotoluene (26DA), and so on.6,7 Several existing experiments showed that the amino metabolites of TNT could accumulate in the tissue with little transformation, and several metabolites are as toxic as TNT or even more toxic.8,9 By exposing catfish to 14C-labeled TNT in water, Ownby et al10 found that TNT only accounted for 4.9% of the total radioactivity after 8h exposure, meanwhile its metabolites could accumulate to a greater extent than the parent compound. In addition, the bioconcentration factors (BCF) of TNT for fishes and aquatic invertebrates are range from 0.31 to 9.71 ml·g-1, lower than those of monoamino metabolites (10.2~13.1 ml·g-1).11 Therefore, the metabolites of TNT should also be taken into account when we try to explore the mechanism of TNT toxicity. The toxicity of TNT and its metabolites in organisms is obviously related to the process of adsorption, distribution, metabolism, and excretion (ADME). So far, we know relatively little about the metabolism process of TNT, although two related mechanisms have been proposed, the oxidative stress and covalent adducts of metabolites binding to proteins and DNA.12,13 The oxidative stress is mainly caused by the formation of nitro anion radicals and their redox cycling, which is initiated by the one-electron reduction and the commonest cause of TNT-induced cataracts.14 Covalent adducts are mainly attributed to the formation of nitroso and/or hydroxylamino metabolites in the process of two-electron reduction, which may modify the structures of proteins and DNA or induce the methemoglobin.15 However, to our knowledge, the study on the interaction of TNT and its metabolites with lipid membranes is rare, in spite of its vital importance to understand their adsorption, distribution, and excretion process in organisms, except for two studies of Golius et al.16,17 They mainly calculated the partition coefficients for several explosives (such as

3



TNT, dinitrotoluene, dinitroanisole, and so on). The detailed information about the interactions between TNT (including its metabolites) and lipid membranes was not provided yet. The lipid membrane is the first barrier for living cells, and the interactions of xenobiotic chemicals with the cell membrane are unavoidable when they enter into organisms, gaining a clear insight into the xenobioticmembrane interactions would help us to understand the ADME and toxic mechanisms of xenobiotics. For example, Qiao and co-workers explained the reduced acute toxicity of functional fullerenes by comparing the translocation of C60 and C60(OH)20 across a model membrane.18

Figure 1. Schematic structure of TNT and its metabolites in this work.

Several experimental techniques, such as ultracentrifugation, solid phase microextraction, equilibrium dialysis, have been utilized to estimate the partition and penetration properties of chemicals and their dependence on concentration, temperature and pressure.19 However, obtaining a detailed molecular level understanding of chemical-membrane interactions is still a challenge for the pure experiment. Computer simulation techniques, such as molecular dynamics (MD) simulations, are powerful tools to gain atomic insights into the interactions of membrane with chemical molecules and relevant permeation process. MD simulations have been successfully employed to study the permeation and distribution of several molecules, as well as how these molecules affect the properties of the lipid bilayer.20-24 Consequently, to deeply understand the ADME process and toxicology mechanism of TNT and its metabolites (Figure 1), in this work, we explored the interactions between them and lipid membranes, and the effect of them on membrane properties, mainly using all-atom MD simulations and umbrella sampling (US) technique.  METHODS

4


Page 4 of 29

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MD simulations All MD simulations in our work were performed by the GROMACS v5.1.2 program suite.25 The bilayer lipid model containing 128 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) molecules were constructed by the CHARMM-GUI,26 and the system was hydrated by around 13,000 TIP3P water molecules. The ionic concentration was set to 0.15 M NaCl. The CHARMM36 force field, which has been proved to be a good choice to compute the permeability of small molecules, was applied to describe all POPC lipids.27,28 Parameters of small molecules (2A, 4A, 24DA, 26DA and TNT) were obtained from the CGenFF force field.29 In particular, the parameters of 4HA, which is not included in standard CGenFF force field, were calculated by the ffTK program.30 We calculated the 1-octanol-water partition coefficient to validate the force field parameters of small molecules. To build the initial models, all small molecules were randomly inserted in the region within 15 Å from the bilayer center along z direction by the Packmol package.31 The final size of simulation box is c.a. 6.7 6.712.2 nm3 and the periodic boundary condition was introduced. To remove the unreasonable contacts, the steepest-descent algorithm was firstly used to minimize energy, and the system was subsequently equilibrated for 1 ns. The NPT ensemble was chosen for all simulations. The systems and corresponding simulation time are summarized in Table S2. For the unbiased MD simulations, the first 80 ns were used for equilibrium unless specified. A velocity rescale thermostat32 was used to maintain temperature at 310 K with a coupling constant of 0.2 ps. Semi-isotropic pressure at 1 bar was controlled by the Parrinello-Rahman barostat33 with a coupling constant of 1 ps. Electrostatic interactions were calculated using the Particle mesh Ewald (PME) algorithm34 with the real space cutoff of 1.2 nm, while the van der Waals interactions were switched to zero smoothly between 1.0 and 1.2 nm. The integration time step was set to 2 fs and the neighbor list was updated after every 10 steps with a 1.2 nm cutoff. The LINCS algorithm35 was used to constrain the bonds involving hydrogen, while the SETTLE algorithm36 was chosen to constrain the bond lengths and angles of water. The trajectory was stored every 20 ps. The snapshots were produced by the VMD 1.9.2.37 PMF calculations Umbrella Sampling (US)38 method was utilized to calculate the potential of mean force (PMF) for the partitioning of molecules considered here into the lipid bilayer. A total of 32 umbrella windows were 5



defined along the bilayer normal ranging from z=0 (bilayer center) to z=4 nm (bulk water). The spacing of windows was set to 0.2 nm in bulk water and 0.1 nm in the bilayer. The center of mass for the molecules was restrained by a harmonic potential with a force constant of 1000 kJ·mol-1·nm-2. Each window was run for 50 ns, totaling 1.6 μs per molecule, and the trajectory was recorded per 20 ps. The PMF profiles were constructed by the weighted histogram analysis (WHAM),39 using the last 20 ns of MD for analysis. In order to obtain an estimation for statistical error for PMF profiles, the Bayesian bootstrap analysis (N=100) was utilized.  RESULTS AND DISCUSSION Single-molecule simulations PMF profiles. In order to obtain the detailed information about the absorption process for TNT and its metabolites, we firstly calculated the Gibbs free energies of penetration from the water phase to the lipid membrane interior at different z-positions along the bilayer normal (i.e., PMF profiles). The corresponding profiles for all species considered are shown in Figure 2. PMF profiles confirm that the membrane absorption of TNT and metabolites is a spontaneous process. There is a distinct drop of free energies when these compounds enter into the bilayer from the water phase, which can be regarded as a driving force (ΔGDF) for the membrane absorption of TNT and metabolites. The size of ΔGDF is closely related to the functional group of molecule, as shown in Figure 2. Specifically, the hydroxylamine metabolite 4HA has the largest driving force for the membrane absorption (ΔGDF = -16.1 kJ·mol-1), which is about 2-fold larger than those of diamino metabolites (-7.8 kJ·mol-1 for 24DA and -6.1 kJ·mol-1 for 26DA, respectively). It is reasonable that the driving force defined by the ΔGDF can be used to describe the membrane affinity of xenobiotics. Thus, based on Figure 2, a clear order for the membrane affinity of TNT and its metabolites can be obtained, that is, 4HA > 4A > 2A > TNT >> 24DA > 26DA. Location and orientation. The location and orientation of solutes inside the membrane are the key parameters of their interaction with membranes which can affect the metabolism, the structure and function of membrane.19,40 The results of PMF profiles predict that the metabolites can penetrate deeper into the membrane than TNT. Moreover, TNT and its metabolites are likely to accumulate at the region with z = 1.29 ~ 1.51 nm, which corresponds to the interface between the head group and tail of lipids. 6


Page 6 of 29

Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. PMF profiles along the bilayer normal (z) for TNT and its metabolites. (a) TNT, (b) 2A, (c) 4A, (d) 4HA, (e) 24DA, and (f) 26DA.

The metabolites with the metabolized groups at 4-position (4A and 4HA) are closer to the tail portion of lipids (Figure 2c and 2d), this result can be rationalized by the orientation difference as discussed below. The orientation is described by the distribution of the tilt angle (θ) for all compounds at zmin, which is defined by the angle between the vector from the carbon at 4-position to the carbon at methyl and the bilayer normal. The results are shown in Figure 3. It is observed that the plane of the benzene ring tends to be parallel to the bilayer normal for these compounds, and the amino group and hydroxylamine group are more likely attracted by the head groups of lipids compared with the nitro group (Figure 3a~3d). For the metabolites with the metabolized groups at 4-position (4A and 4HA), the nonpolar methyl tends to point to the hydrophobic lipid tails (Figure 3c, 3d), which is more favorable for solutes to enter into the hydrophobic regions in the membrane compared with the 2-position metabolites (2A). This result well corresponds to the observations about the location mentioned above. Specifically, in the case of 24DA (Figure 3e), the two adjacent amino groups compete and both closely interact with head groups, resulting in the most probable

7



tilt angle of about 120°, the similar results are also observed for TNT which has two adjacent nitro groups. It should be stressed that the molecules with similar functional groups (e.g. 2A and 4A) have distinct orientation as well as location, which may be one of the reasons for their different behaviors at high concentration such as aggregation.

Figure 3. The orientation and tilt angle distribution at the position with free energy minimum (zmin). (a) TNT, (b) 2A, (c) 4A, (d) 4HA, (e) 24DA, and (f) 26DA.

Membrane affinity and Penetration properties. Generally, the partition coefficient between water and lipid membrane is also chosen to describe the membrane affinity of xenobiotics.16 Then, this parameter can provide an important theoretical reference to understand the differences of the bioaccumulation of TNT and its metabolites, especially in membranes. Previously, the octanol-water partition coefficient (KOW) was usually used to estimate the bioaccumulation because of the large availability of experimental data and simple measurement for KOW.41 However, since the octanol cannot possess all complexity of a membrane structure, the KOW can only characterizes the hydrophobicity rather than lipophilicity, further not describe the partition of chemicals between the membrane and water bulk accurately.17 A new liposome-water partition coefficient (KLW) was recently suggested to provide an accurate estimation for the membrane 8


Page 8 of 29

Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


affinity of chemicals.41 In MD simulations, the value of KLW can be constructed from PMF profiles according to the following equation,42 ∑𝑛

𝐾LW =

𝑥w(𝑧𝑖) 𝑉 𝑧𝑖 exp ( ―𝛽∆𝐺(𝑧𝑖)) ― 𝑉(𝑧𝑛)exp ( ―𝛽∆𝐺(𝑧𝑛))𝑥 (𝑧 ) w 𝑛

[()

𝑖=1

]

(1)

1

𝑉(𝑧𝑛)exp ( ―𝛽∆𝐺(𝑧𝑛))𝑥 (𝑧 ) w 𝑛

where the system is split to layer 1 to n along z direction and the layer n belong to the water phase. β = 1/ 𝑘𝐵𝑇; kB is the Boltzmann constant; T is temperature (K); 𝑥w(𝑧𝑖) =

𝑛w(𝑧𝑖) 𝑛tot w

; 𝑛w(𝑧𝑖) is the number of water

molecules in layer i and 𝑛tot is the total number of water molecules in the system. It is assumed that a w fraction of the free energy in each layer comes from the water content in this layer, and this contribution is subtracted from the free energy. Here all the water molecules in the system are taken into account. This equation considers the whole free energy profile and avoids the definition of limit between water phase and lipid phase. The calculated results of log KLW and the corresponding log KOW are listed in Table 1, the data of log KOW were obtained from previous work except for 4HA which is from our simulations. It is presented that the values of log KLW is very close to that of log KOW, the average difference between log KLW and log KOW is only 0.12 and the largest difference (0.31 log unit) is obtained from4A. Previous study41 into plenty of neutral organic chemicals showed that the values of log KLW and log KOW exhibit a 1:1 agreement with the average difference of 0.4 log unit. Our result is in good agreement with this conclusion, which shows that our simulations give good estimation of the log KLW and also provide evidence for the validation of our MD parameters and models. It is observed that the order of log KLW is not exactly the same as the order of log KOW, e.g., log KLW (4A) > log KLW (2A) while log KOW (4A) ≈ log KOW (2A), and we believe the order of log KLW can better describe the membrane affinity of TNT and its metabolites. Based on the results of log KLW in Table 1, the order of their membrane affinity should be 4HA > 4A > 2A > TNT >> 24DA > 26DA, which is in complete agreement with the result obtained from the ΔGDF. Therefore, both parameters (ΔGDF and KLW) can give an acceptable estimation of the membrane affinity for xenobiotics. In general, the adsorption process of xenobiotics in membrane is rationalized by a free-energy profile (i.e., PMF profile) along the lipid bilayer normal, which can be reliably reconstructed from biased MD simulations. Here, all PMF profiles have the similar shape, the characteristic of which is an energetic minimum (zmin) in the curve, as shown in Figure 2. Based on Eq. (1), we can see that the KLW is decided by the whole PMF 9



curve, but the free energy with smaller values contribute more to the KLW, and the minimum of the free energy ΔGmin (i.e., ΔGDF in our simulation) makes the largest contribution to the KLW. That’s why the order of ΔGDF can be used to quantitatively estimate the order of membrane affinity in our simulation. We next estimate bioaccumulation of these compounds by calculating the bioconcentration factor (BCF) based on the values of KOW and KLW. BCF is the ratio of a chemical concentration in the organism to that in the surrounding water, and it is suggested that a directly proportional relationship exists between BCF and octanol-water partition coefficient (KOW).43 In present study, the following equation proposed by Meylan et.al.44 is used to predict the BCF: log BCF = 0.86log KOW ―0.39

(2)

we assume that this equation is also valid for KLW, the values of BCF based on KOW and KLW (named as BCFOW and BCFLW respectively) are collected in Table1 and compared with the experimental results. It is shown that the predicted values are generally larger than experimental values. The results of BCFOW and BCFLW do not show significant difference except 4A, it is not surprising because KLW and KOW have similar values for these molecules. In the case of 4A, BCFLW (33.1) is obviously greater than BCFOW (17.9), leading to the different order of predicted BCFs, e.g., BCFLW (4A) > BCFLW (2A) while BCFOW (4A) < BCFOW (2A). In previous experiment, Conder et.al.9 found that the BCFs for 4A (12.41) was slightly larger than 2A (10.22), in qualitative agreement with the results of BCFLW. The results suggest that the liposome-water partition coefficient can give a reasonable estimation on the bioaccumulation of xenobiotics. The distribution and excretion in ADME processes of xenobiotics mainly depend on the penetration of xenobiotics and the resulting metabolites through cell membranes. In such processes, the parameter of KLW only describes the membrane affinity of a substance. Therefore, to estimate the capability of TNT and its metabolites penetration across the overall membrane and the time needed for this process, we further calculated the permeability coefficient (P) and the transmembrane time (τ). According to the inhomogeneous solubility-diffusion model,45 the coefficient P can be calculated from the integral of the ratio between exp(ΔG(z)/RT) and D(z), as shown in Eq.(3),

[

4nm 𝑒

P = ∫ ―4nm

Δ𝐺(𝑧) 𝑅𝑇

𝐷(𝑧)

𝑑𝑧

]

―1

(3)

10


Page 10 of 29

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where R is the gas constant and T is the temperature (K), D(z) is the diffusion coefficients at position z along bilayer normal. The integration extremes -4 nm and 4 nm correspond to two sides of the bilayer membrane. The calculated P of TNT and its metabolites are listed in Table 1. It is showed that the order of the P for TNT and its metabolites is 2A > 4HA > 4A >> 24DA > 26DA > TNT. It is observed that the metabolites have larger permeability than TNT, implying their greater ability to cross membranes. Interestingly, the order of P is exactly the same as the reverse order of pure barrier (ΔGPB, equal to the free energy difference from water phase to the position of free energy maximum) in PMF profiles, for example, TNT with the smallest P (0.017 cm·s-1) has the highest pure barrier (ΔGPB = 15.9 kJ·mol-1, as shown in Figure 2a). This result can also be rationalized from Eq. (3), in which the ΔG with the positive value (i.e., ΔGPB) makes the largest contribution to the permeability coefficient, considering the characteristics of exponential function. In addition, the average transmembrane time for xenobiotics was also estimated by the assumption of Brownian motion as18,46 1

4nm

𝑦

τ = 𝐷𝑧∫ ―4nm𝑒∆𝐺(𝑦)/𝑅𝑇𝑑𝑦∫ ―4nm𝑒 ―∆𝐺(𝑧)/𝑅𝑇𝑑𝑧

(4)

Where Dz is the diffusion coefficient in the transmembrane direction. As equation (4) is a combination of homogeneous and inhomogeneous diffusion model, with an effective Dz instead of a position dependent Dz(z), it only serve for an order-of-magnitude calculation. Dz is obtained by averaging the value of Dz(z) over the region z = 1.3-1.5 nm, where the molecules most likely inhabit during the translocation process. The results are also presented in Table 1. The final order of the τ for TNT and its metabolites is TNT > 4A >> 4HA > 24DA > 2A > 26DA. The transmembrane time of TNT is generally an order of magnitude larger than its metabolites, indicating that TNT is more likely to be slowed by the membrane, which agrees with the result of permeability coefficients. In addition, from the results of τ in Table 1, we can know that 4A (τ = 110 s) is more likely to be slowed by the membrane than 2A (τ = 23 s). We can explain this phenomenon by PMF profiles. Based on Eq.(4), the τ relies on the total free energy barrier in the membrane center (ΔGPB + ΔGDF). In detail, the total barrier for 4A (ΔGPB + ΔGDF = 27.6 kJ·mol-1) is 4.3 kJ·mol-1 higher than that of 2A (22.7 kJ·mol-1), indicating that 4A needs more time to cross the higher barrier if it would penetrate the lipid bilayer.

11



Page 12 of 29

Table 1. Summary of the results calculated from single-molecule simulations and related literatures. ΔGDF (driving force for the membrane absorption, kJ/mol), log KLW (liposome-water partition coefficient), log KOW (octanol-water partition coefficient), BCFLW (The predicted BCF based on log KLW, ml g-1), P (permeability coefficient, cm/s), τ (transmembrane time, s), BCFOW (The predicted BCF based on log KOW, ml g-1). Results from literatures

Results calculated from our simulations Compound

a.

ΔGDF

log KLW

BCFLWb

P

τ

log KOW

BCFOWa

BCF(expt)c

TNT

-12.0

1.61 (0.26)

9.9

0.017

430

1.6c

9.7

0.31-9.71

2A

-14.7

2.01 (0.11)

21.8

0.99

23

1.94c

19.0

10.22

4A

-15.9

2.22 (0.11)

33.1

0.24

110

1.91c

17.9

12.41

24DA

-6.8

0.79 (0.27)

1.95

0.069

24

0.7c

1.63

2.75

26DA

-5.5

0.54 (0.24)

1.19

0.044

20

0.67c

1.54

-

4HA

-16.1

2.49 (0.26)

56.4

0.67

62

2.6d

70.1

-

Predicted by the equation log BCF=0.86 log KOW-0.39 from Meylan et al.44

b. Using the same equation as a, and log KOW is replaced by log KLW. c.

From ref 11.

d. Calculated based on our force field parameters.

The results of single-molecule simulations give insight into the process of TNT and its metabolites to cross the cell membranes. We first obtain the PMF profiles from US calculation, based on this curve a basic understanding of the transmembrane process is provided and several valuable parameters are calculated. The liposome-water partition coefficients give the membrane affinity of these compounds, estimating their capability to enter into membrane from water phase. From the permeability coefficient and transmembrane time, the ability of TNT and its metabolites to transfer across the overall membrane and enter into cells is quantitatively assessed. In addition, their behaviors inside membrane are also explored by obtaining the orientation and distribution. These data are helpful for us to understand the overall transmembrane process of TNT and its metabolites and the interactions of them with cell membranes in this significant process. Multi-molecule simulations Membrane properties. In the above study, we obtained a comparison of the bioaccumulation ability of TNT and its metabolites in lipid bilayers by US calculations of single molecules. The results present that 12


Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TNT and its metabolites can enter into membrane spontaneously and accumulate inside the bilayers, which can further lead to the change of membrane properties. Next, we will focus on the effect of multi-molecules accumulation (i.e., concentration effects) on membrane properties and its causes for TNT and its reduced metabolites, because it may be related to their toxic mechanisms. It should be noted that the high concentration (may be much higher than the real physiological condition) is used in our work to distinguish the concentration effect from statistical noise or thermodynamically fluctuation.24 Moreover, Kopec et al47 pointed out that the concentration in MD simulations could not be directly related to the physiological concentration because the simulation can only consider a nanometer-sized membrane patch, which might exist in a local region or at a specific time.

Figure 4. The structural properties of membrane at different concentrations of TNT and its metabolites. (a) Thickness, (b) Area per lipid, (c) Average order parameter, and (d) Lateral diffusion coefficient. The continuous connecting curves do not represent physical relationships, and are just guides to the eye.

Here, three different numbers (20, 40 and 60) of target molecules are chosen to check the concentration effects on POPC membranes, and the corresponding membrane properties are calculated and compared in Figure 4 (The detailed data are given in Table S3). Firstly, for the thickness of lipid bilayers (TLB), it is 13



Page 14 of 29

obtained as the average distance between the phosphate atoms of two monolayers. The results in Figure 4a show that all TLBs have a slow decrease with the increase of concentrations of target compounds. The largest decline, from 3.89 nm to 3.68 nm for 60 4HA in membranes, is also less than 5.4%. For the low concentrations (i.e., 20 molecules in membranes), the largest drop of TLB is about 1.2%. Secondly, for the area per lipid (APL), it is calculated by the following equation, 𝐿𝑥𝐿𝑦

APL = 2 𝑁𝑙𝑖𝑝

(5)

where Lx and Ly are the box length in X and Y direction, respectively, and Nlip is total number of lipids in the bilayer. As shown in Figure 4b, the APL generally increases with the concentration of TNT and its metabolites. The maximum increase appears in the system with 60 4HA molecules, in which the APL increases by about 14%, from 64.7 Å2 to 73.8 Å2. Both the augmentation of APL and decline of TLB suggest that the membranes should become more disordered when TNT and its metabolites interact with hydrocarbon tails. Therefore, we estimated the order degree of different systems by the deuterium order parameter (SCD) of hydrocarbon tails, which is defined as: 1

𝑆CD = 2〈3 cos2𝜃 ― 1〉

(6)

where  is the angle between the C-D (C-H in the present simulations) vector and the bilayer normal. The brackets indicate the ensemble average. The order parameter SCD for sn-1 chains in the lipid tails are given in Figure S1. For the pure bilayers, our calculated order parameters are in agreement with previous ∑𝑆𝐶𝐷

simulation or experimental studies.48,49 The average order parameters 〈𝑆CD〉 over all C-H vectors (

𝑛

) for

all species are drawn in Figure 4c. The values of 〈𝑆𝐶𝐷〉 for 4A, 4HA, 24DA and 26DA show an obvious decrease with the increase of their concentrations. The biggest variation is found for the systems with 60 molecules. Interestingly, for TNT and 2A, values of 〈𝑆𝐶𝐷〉 do not change significantly when the number of TNT and 2A in the membrane exceeds 40. 〈𝑆𝐶𝐷〉 of 60 2A system is even larger than that of the 40 2A system, which indicates that the degree of ordering of lipid tails slightly increases. Moreover, the degree of ordering should also be influenced by the distribution of compounds in membranes. In the next section, it is suggested that TNT and 2A actually forms the aggregations in the membrane which influence their distributions. Thus, we can know that the interactions of TNT and its metabolites with the lipids reduce the 14


Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


orderliness of lipid tails, which, in turn, results in the increase of APL, and this process would also be affected by the distribution of species. Finally, to check the lateral dynamic of membranes, we calculate the lateral diffusion coefficient (LDC) of the phosphorus atom in the lipid head group from the slope of meansquare displacements (MSD) in the XY-plane (Figure 4d). The LDC value of pure POPC membrane is 0.9110-7 cm2·s-1 in our simulation, which is in good agreement with previous experimental result.50 The results in Figure 4d suggest that the LDCs are not significantly affected by the penetration of TNT and its metabolites, in spite of the high concentration which has been used. It should be noted that the concentration of target chemicals used in our simulations may be much higher than that of physiological condition, so the change of membrane properties in the actual situation should be much smaller. Whereas, the high concentration used in our simulation is helpful to observe the variation tendency with the concentration. Moreover, such a high concentration may also be achieved in some cases. Lachance et al.8 found that the highest concentration of 2A in earthworm tissues is 3.0 μmol·g-1, corresponds to about 3.0 mM on assumption of the density of 1g·mL-1. The highest concentration of 2A in our simulations is ~100 mM, quantitatively comparable with the physiological concentration in earthworm tissues, if we consider that 2A more likely accumulate in lipid membrane than other tissues. Distributions in membranes. The previous analysis pointed out that the concentration effect of TNT and metabolites also depends on their distributions in membranes. Thus, we provide the final snapshots of all systems in Figure 5 and S2, Figure 5 only shows two representative examples, that is, the monoamino metabolite 2A and 4A. For all systems, the species prefer to randomly locate below the head group, well corresponding to the results of PMF profiles. Interestingly, for the system with 60 molecules of 2A, a membrane-spanning aggregation is formed; however, for the similar 4A system, the aggregation is not observed, as shown in Figure 5. We extend the simulations to 240 ns and the result does not change. For the systems with 40 and 60 TNT molecules, a similar aggregation is also found, respectively. However, we do not observe any membrane-spanning aggregations for other TNT metabolites, even though the identical high concentrations are used (Figure S2). Our other simulation study explored the origin of TNT aggregation (under review), the results show that the VDW interaction between nitro groups of TNT molecules promotes the formation of aggregations. We also found that the existing TNT molecules inside

15



the membrane will decrease the central barrier for subsequent TNT to cross the membrane center, allowing the aggregation to span the overall membrane. To quantitatively describe the distributions of TNT and metabolites, the mass density in membranes are calculated and shown in Figure 6. It is observed that the species tend to occupy the region between the head group and the hydrocarbon tail, which is in good agreement with the results of PMF profiles and the snapshots in Figure 5 and S2. In addition, the peak of mass density generally increases with the increase of species concentrations, except for 2A in Figure 6b. It is presented that the peak density does not increase when the number of 2A increased from 40 to 60. However, a considerable number of 2A molecules remain in the central region of the membrane, that is, between the two peaks. As shown in Figure 5c, those molecules form a membrane-spanning aggregation. The same phenomena is also obtained for TNT molecules (Figure 6a). A little difference is that when the number of TNT molecules exceeds 40, they accumulate in the region between the two peaks. A detailed comparison shows that when the number of species increases from 20 to 40, the peak density of TNT has the minimum increase, compared with its metabolites. It can be attributed to the formation of membrane-spanning aggregations.

Figure 5. Final snapshots of (a) 20 2A, (b) 40 2A, (c) 60 4A, (d) 20 4A, (e) 40 4A, and (f) 60 4A system. The lipid tails are omitted for clarify. The color rules are: lipid heads (yellow), water (cyan), 2A (pink), 4A (magenta). The periodicity is removed to maintain the integrity of molecules.

Finally, from Figure 6e and 6f, it can be seen that, compared with the other species, the mass density curves for the diamino metabolites 24DA and 26DA have a slight difference, that is, beyond the region -2 nm to 2 16


Page 16 of 29

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nm (corresponding to the lipid bilayer region), their concentration is not zero. It is suggested that the diamino metabolites of TNT readily enter into the water phase, even though they are initially placed inside the membrane, which is also supported by Figure S2g~l. Based on the above calculated liposome-water partition coefficients (KLW), 24DA and 26DA have the minimum membrane affinity. Moreover, their PMF profiles also show that they possess a very low free energy barrier (7.8 kJ·mol-1 and 6.1 kJ·mol-1 for 24DA and 26DA, respectively) from the membrane interior to the water phase. These results fully demonstrate that the diamino metabolites of TNT are not easily accumulated in cell membranes.

Figure 6. Mass density profiles for TNT and its metabolites at different concentrations. (a) TNT, (b) 2A, (c) 4A, (d) 4HA, (e) 24DA, (f) 26DA.

Water penetration. The structure of membranes is actually influenced by TNT and its metabolites. In general, the change of membrane structure will affect its function. Due to the amphiphilic character of lipid molecules, the membrane acts as a barrier of separating the cytosol from extracellular environment.

17



Specifically, the membrane interior is hydrophobic and will prevent water molecules from staying in this region. Therefore we check the permeability of water in the POPC membrane with the different numbers of TNT and metabolites. The main method is to estimate the number of water molecules entering into the membrane interior. With the center of the bilayer as the origin (z = 0 nm), the number of water molecules penetrating into three membrane regions was calculated, that is, Region I of 0.9 nm < |z| < 1.3 nm, Region II of 0.5 nm < |z| < 0.9 nm and Region III of |z| < 0.5 nm. According to the regions defined by Marrink and Berendsen,45 Region I and II correspond the high tail density region and Region III corresponds to the low tail density region. The results of water numbers counted in three regions are given in Figure 7. The vast majority of water molecules can generally penetrate the region of 0.9 nm ~ 1.3 nm (or -0.9 nm ~ 1.3 nm). In the Region II of 0.5 nm < |z| < 0.9 nm, very few water molecules have been found, except for the systems with the high TNT accumulation. As shown in Figure 7a, the number of water molecules entering the membrane shows an increase trend with the increasing adsorptions of TNT and its metabolites. Khondker et al.51 also found the similar trend in the system of caffeine molecules interacting membranes, which would bring the loss of membrane fluidity. For the most hydrophobic Region III (closest to the center of membranes), only the system with 60 TNT molecules finds traces of water molecules (Figure 7c). In the middle Region II, it is also found that the systems with high concentrations of TNT (40 and 60 TNTs) could induce more water molecules to enter the membrane (Figure 7b). In fact, when the number of TNT molecules exceeds 40, an aggregation is formed in the lipid bilayer (Figure S2e and S2f), and it can lead to the formation of water flow through the overall lipid bilayer along the TNT aggregate (data not shown). In addition, the monoamino metabolite of TNT (2A) also forms an aggregation in the membrane when the number of 2A molecules reaches 60. However, the ability of 2A aggregates to induce water molecules into membranes is much weaker (Figure 7).

18


Page 18 of 29

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. (a) The number of water molecules in the different defined regions of the lipid bilayer interior vs. different concentrations. The region is defined as follow: Region I (0.9 nm< |z|

Molecular Simulation Studies on the Interactions of 2,4,6

Recommend Documents