A comprehensive model for chemical bioavailability and toxicity of organic chemicals based on first principles

Introduction

There is a large number of chemicals, both natural and synthetic that can bind to the aryl hydrocarbon receptor (AhR), sometimes referred to as the “dioxin receptor” (Giesy et al., 1994a,b). Most research into the toxicity of AhR-agonists and other chemicals has focused on assays with actual biomaterial, which can be time-consuming and require use of live animals. Alternatively, modern molecular techniques and in silico studies of quantitative structure-activity relationship (QSAR) methods are being used more and more to predict biological activities of organic molecules (Safe, 1993; Yang et al., 2009, 2010). Such approaches are, however, limited in their ability to explain the inherent nature of bioactivity (viz., potential toxicity) of chemicals because of their case- or compound-specific masking effects (Safe, 1993; Chana et al., 2002). This method is not always accurate and uses information on the physical structure of molecules to calculate toxicity.

Polychlorinated biphenyls (PCBs) were manufactured and used as electrical insulating liquids from 1929 to the late 1970s, when they were voluntarily withdrawn by some manufacturers and subsequently banned in Europe and North America, because of their environmentally persistent, bioaccumulative, and potential toxic effects on wildlife and humans (DiGiovanni et al., 1987). Some congeners of PCBs, which comprise a group of 209 possible compounds with two chlorinated 6-member carbon (C) rings connected by a C-C bond (Giesy and Kannan, 1998) are potent, chronic toxins (Giesy et al., 2006). Previous researches have shown that, in animals or in vitro, PCBs are potentially carcinogenic (DiGiovanni et al., 1987), neurotoxic (Seegal, 1996), and also affect the endocrine system (Birnbaum, 1994).

Toxicity of PCBs is known to be dependent on the pattern of chlorine substitution, viz. their numbers and positions allocated, on the phenyl rings (Giesy et al., 2006). For example, a subgroup of 12 PCBs (non-ortho- or mono-ortho-chlorinated) exhibit “dioxin-like” activities associated with the AhR. It has also been shown that the energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is an indication of stability in similar molecules (Lynam et al., 1998) where a larger gap implies a more stable molecule with respect to reactions with biomolecules. The relative stability of more chlorinated PCBs with greater lipophilicity results in widespread distribution in various environments followed by long-term accumulation in wildlife and humans (Alcock et al., 1998). Despite these extensive previous studies, toxicity of PCBs is known with limited accuracy and this makes PCBs an ideal candidate to test the proposed model for the toxicity and bioaccumulation behaviors of AhR-agonists and antagonists (Newsted and Giesy, 1987).

Here, an approach to predict toxicities of AhR-agonists from their physicochemical nature, namely their electronic characteristics, is presented. PCBs were selected as a model system to prove the concept of explaining toxicity based on chemical and electronic structure. The proposed first principles model allows for a fundamental understanding of toxic potency of individual PCB congeners but can be applied to other AhR-agonists, including both natural and synthetic compounds, encompassing both chlorinated and brominated or chloro-bromo-compounds. There are many natural and synthetic chemicals that can interact with the AhR (Giesy et al., 1994a). This model will also allow for the rapid screening and classification of tens of thousands prospective AhR-active chemicals without the need to use animal testing and thus help industry develop the chemical products with minimal environmental/toxic impact while retaining useful properties.

Materials and Methods

Density Functional Theory Calculations

Ab initio DFT calculations were performed by use of Gaussian 03 software (Frisch et al., 2004). All calculations employed the 3-21G basis set. The exchange correlation functional and basis sets used were B3LYP and 6-21G. After geometry optimizations, density of states (DOS) and electric DM of each PCB congener were determined by calculations of occupied and virtual (unoccupied) molecular orbital (MO) energies. Mulliken charge population analysis and isosurface visualizations were also performed.

Chemicals and Sample Preparations

For X-ray spectroscopic measurements, selected PCB congeners including PCB-101, PCB-105, PCB-118, and PCB-153 were obtained from Wellington Laboratories (Guelph, Canada). One set of these congeners was then exposed to a rainbow trout microsome and incubated at 37°C for 3 h. The microsome was then extracted, leaving only metabolites of the reaction. Those metabolites were then dissolved in 200 μL hexane. Films of this hexane solution were fabricated on Si (100) wafers using spin-coating at 800 rpm for 20 s.

NEXAFS and XES Measurements

The unoccupied electronic structure of PCB molecules before and after exposure to microsomes was investigated by employing Near Edge X-Ray Absorption Fine Structure (NEXAFS) spectroscopy. Measurements of C 1s NEXAFS were made by use of the Spherical Grating Monochromator (SGM) beam line (Regier et al., 2007) at the Canadian Light Source of the University of Saskatchewan (Saskatoon, SK, Canada). Spectra were measured in total electron yield (TEY) mode. The C 1s NEXAFS spectra were energy calibrated with the π* (C=C) transition at 285.5 eV of highly ordered pyrolytic graphite (HOPG). Samples were tested for radiation damage because of the sensitivity of carbon-based materials to X-ray radiation. This was performed by comparing results of five repeated measurements of the desired energy range for 1–2 min on the same sample spot. To reduce radiation damage, concurrent measurements were taken at different locations of the same sample. Measurements made using NEXAFS spectra were intensity-normalized to the incoming photon flux as recorded by a photodiode. Non-resonant C Kα (2p → 1 s transition) XES were taken at Beam line 8.0.1 of the Advanced Light Source in Berkeley, CA (Jia et al., 1995), with excitation energy of 310 eV for all samples.

Results and Discussion

First Principles Model of Toxicity

PCBs consist of two benzene rings bonded at opposite carbon sites, with varying number and location of chlorine termination. The molecular structure of PCBs and the chlorine termination, indicating nomenclature for the congeners, are shown in Figure 1. The 12 dioxin-like PCBs out of a total of 209 congeners have a few structural similarities. These similarity factors are: (i) co-planarity between two benzene rings, (ii) either zero or one Cl-substituted ortho position, (iii) both para positions chlorinated, and (iv) two or more meta positions chlorinated.

Figure 1. HOMO-LUMO gap of PCB congeners. Calculated HOMO-LUMO gap of all PCB congeners differentiated by number of ortho-chlorination sites. Below threshold of 4.87 eV are the 12 dioxin-like PCBs, above the threshold are non-dioxins. Note that coplanar PCBs have (on average) lesser band-gaps. PCB-74 (denote as A) is considered a possible dioxin, due to its band-gap below the threshold. Torsional angle, α is defined as “twist” around center axis, as depicted in inset.

The dioxin-like PCBs are known to activate the AhR, otherwise referred as the “dioxin receptor.” Since this involves chemical bonding through sharing or transferring of electrons, the above similarities must be associated with the electronic configuration of PCBs (Farmahin et al., 2013). Especially, HOMO and LUMO levels determine how a chemical molecule shares its valence electron in the occupied molecular orbitals and donates (or accepts) electrons to ligand, and thus the HOMO-LUMO energy gap has been suggested as an indicator of the stability of molecules (Lynam et al., 1998). Density functional theory (DFT) numerical calculations were made to obtain HOMO-LUMO energy gaps for all 209 PCB congeners. Calculated HOMO-LUMO gaps, calculated as a function of torsional angle between two benzene rings (α) are given in Figure 1.

Most coplanar PCBs (α = 0°) had lesser HOMO-LUMO gap energies than did non-coplanar congeners. The 12 dioxin-like PCB congeners are among those having lesser energy gaps. PCB-77 (3,3′,4,4′-tetrachlorobiphenyl) has an energy gap of 4.87 eV, which is proposed as a threshold energy gap for toxic potency of dioxin-like PCBs (see Figure 1). The results of these predictions are consistent with those of other studies showing that coplanar PCBs are more toxic through AhR-mediated mechanisms than are the di-ortho-substituted PCBs (Safe, 1993; Van den Berg et al., 2006). There was a relatively narrow distribution of HOMO-LUMO gaps among coplanar PCBs, with a range of 0.35 eV for all coplanar and only 0.09 eV for dioxin-like congeners. The range was wider for the other non-coplanar PCBs (1.2 eV). A larger range of energy-gaps for non-coplanar PCBs with various α angles and no PCB with a small α angle near the coplanar geometry suggest that dioxin-like PCBs are in a metastable state. Therefore, it is expected that dioxin-like PCBs would transform to non-coplanar configurations with larger HOMO-LUMO gaps and stabilize themselves by changing the configurations of chlorine atoms during reaction with the ligand binding domain (LBD) of the AhR protein.

The HOMO-LUMO gap for PCB-74 (2,4,4′,5-tetrachlorobiphenyl) is less than the suggested threshold for significant binding affinity to the LBD of the AhR. PCB-74 is not included in a subgroup of 12 dioxin-like PCBs to be classified as having two or more meta-chlorine atoms but it does satisfy the other similarity factors. The model developed here includes PCB-74 as a dioxin-like PCB, but is classified as a “possible” dioxin-like PCB.

Based on DFT calculations, the HOMO-LUMO energy gap is an overarching indicator of potency of AhR-mediated effects of PCBs that can distinguish the dioxin-like and non-dioxin-like molecules (Similarity 1). The trend observed in HOMO-LUMO gap and potency of PCB congeners is not consistent with the values for 2,3,7,8-TCDD toxic equivalency factors (TEFs) of PCBs (Table 1). This means that a small HOMO-LUMO energy gap that is less than the threshold is a necessary, but not a sufficient condition for significant potency of molecules and there must be other factors to consider when relating toxicity within the dioxin-like subgroup.

Table 1. Structural and electronic configuration parameters of the dioxin-like PCBs and the toxicity ranking predicted by first principles model.

HOMO and LUMO isosurfaces of the most potent dioxin-like PCB-126 (3,3′,4,4′,5-pentachlorobiphenyl) and those for non dioxin-like PCB-107 (2,3,3′,4′,5-pentachlorobiphenyl) did not seem to demonstrate any AhR-mediated potency (Figure 2). Although there were no significant differences in spatial distribution of HOMO and LUMO states between two molecules with different chlorination configuration, the chlorinated sites contributed differently to the HOMO and LUMO states. For this reason another site-specific factor representing how sensitive each chlorinated site is to electron transfer (accepting or donating an electron) when a PCB molecule is interacting with the AhR was considered for inclusion in the predictive model.

Figure 2. HOMO and LUMO isosurfaces of PCBs. The HOMO and LUMO isosurfaces of PCB 126 and PCB-107.

Examination of the Cl-terminated sites of the 12 dioxin-like PCBs, Δ₊(Δ₋)is the absolute difference in partial charge at the Cl-terminated sites between neutral and positively (negatively) ionized molecules. This factor, expressed in units of electron charge, e, is a measure of which sites are most chemically active, and is thus referred to as the “site-specific reactivity” factor. A Mulliken charge population analysis (Frisch et al., 2004) was performed for all dioxin-like PCB congeners in negatively and positively ionized states, as well as an electrically neutral state (Figure 3). Para positions (4 and 4′ sites) were most likely to either accept or donate more of its partial charge. Especially, when positively ionized, para sites have about 50% larger Δ₊ than other chlorinated sites. This observation is consistent with the para positions being most favorable sites for interaction with the AhR, followed by meta sites then ortho sites. Results of the Mulliken analysis were consistent with similarity factors 3 (both para positions chlorinated) and 4 (two or more meta positions chlorinated). Furthermore, Δ₊ was found to be consistently larger than Δ₋ across all sites and PCB molecules, which suggests that for a given interaction PCBs preferentially donate electron charge and become positively ionized.

Figure 3. Average site-specific reactivity of PCBs. Average site-specific reactivity for (A) giving up electron (B) taking electron for the 12 dioxin-like PCBs. Note the para (4 and 4′) sites are most active, then meta (3, 3′, 5, 5′) then ortho (2, 2′, 6, 6′). Also indicated are the average activities of non-ortho and mono-ortho dioxin-like PCBs. In both cases the non-ortho PCBs are more active than mono-ortho.

The partial charge differences (Δ₊ and Δ₋) averaged over all chlorinated sites of non-ortho-substituted PCBs (red dashed lines) are greater than those of mono-ortho-substituted molecules (blue dashed lines). Overall chemical reaction behavior of the dioxin-like PCBs is therefore suppressed when the ortho position is chlorinated. The ortho chlorinated site has more delocalized HOMO states than do meta sites in the same benzene ring (3 and 5 sites). This means that the more the ortho positions are chlorinated, the more difficult it is for a PCB molecule to retain its coplanar structure. This is due to the Coulomb repulsion of delocalized π (occupied) orbitals of the two benzene rings. By considering these effects of the ortho-chlorination, the 12 dioxin-like PCBs could be separated into two subgroups, the more potent non-ortho PCBs and the less potent mono-ortho PCBs (Table 1).

The relative toxic potency of each AhR-active PCB in each subgroup was then calculated. The first principles calculations suggested that binding of dioxin-like PCB with a receptor molecule, such as the AhR is not a random event but is strongly influenced by electronic configurations.

The dipole moment (DM) induced by nonhomogeneous charge distribution in a molecule can be a useful parameter for prediction of toxic potency, but the DM of the 12 dioxin-like PCBs varies widely from 0 to 3.3 D with no clear trend (Chana et al., 2002). Therefore, direction as well as magnitude of the DM should be considered relative to affinity of binding and thus potency of toxicity. It was hypothesized that the toxic potency would be greater when the DM direction with respect to the center axis of the molecule (θ as depicted in Figure 4) maximizes exposure of chemically reactive chlorinated sites to the AhR LBD. The 4 non-ortho PCBs in the more toxic subgroup, PCB-126 (3,3′,4,4′,5-pentachlorobiphenyl), PCB-169 (3,3′4,4′,5,5′-hexachlorobiphenyl), PCB-81 (3,4,4′,5-tetrachlorobiphenyl), and PCB-77 (3,3′,4,4′-tetrachlorobiphenyl) were considered initially. As given in Table 1, PCB-81 has a DM directed along the center axis of a PCB (θ = 0°) and, thus the molecule can align such that the most reactive para sites are preferentially exposed to a neighboring receptor. In the case of PCB-77, the direction of the DM is perpendicular to the center axis (θ = 90°), which results in preferential exposure of the 3 and 3′ meta sites. This explains why PCB-81 has greater toxic potency than does PCB-77. PCB-126 has a DM directed θ ~ 61° below the molecule's center (see Figure 4). Although the DM alignment in this case has lesser site-specific reactivity contributed by para sites, the contribution from more meta chlorinated sites in PCB-126 than in PCB-81 and PCB-77 could compensate for the angle of the DM, producing a greater toxicity in PCB-126.

Figure 4. Direction of dipole moment. PCB-126, showing direction of dipole moment (DM), relative to center axis of molecule (θ ). θ is defined as clockwise the below center axis, 0° < θ < 360°. The DM in this example preferentially allows more active para and meta sites to attach to bio-matter.

Based on this analysis, a function that describes the directional reactivity factor (DRF), analogous to TEF, was derived (Equation 1).

This function represents the difference in partial charge, or site-specific reactivity (Δ₊), for each of the para and meta sites weighted by the direction and magnitude of the DM which maximize or minimize exposure of chlorine sites. Δ_+n is the normalized site-specific reactivity for a para or meta site n when removing one electron from the molecule. All six sites were summed regardless of chlorination where terms from chlorinated sites are positive (contribute to reactivity due to the electronegative Cl) and those from hydrogen sites are negative.

The site-specific reactivity was normalized to a value of 1.0 within the non-ortho and mono-ortho subgroups. $\left|{\stackrel{\rightharpoonup }{\text{p}}}_{\text{n}}\right|$ is the magnitude of the DM which is also normalized to 1 within the subgroups, such that it is of equivalent scale to reactivity. cosϕ_n represents the component of DM along the direction of the nth site's reactivity which is taken to act along the direction of the Cl-C (or H-C) bond for each para or meta site. Since all dioxin-like PCBs have an identical phenyl ring structure, the angular components ϕ_n are constant for all non-ortho or mono-ortho molecules (see Figure 4). The absolute value of the cosine was used because the DM acts as a weight for the reactivity only and should not produce negative values. That is, chlorinated sites always contribute positively to the DRF. The possible reduction of chemical reactivity due to dipole direction was taken into account by including hydrogen para and meta sites with normalized negative contribution.

This model, based on first-principles, was determined to work for non-ortho PCBs and mono-ortho-substituted, or dioxin-like PCBs (Table 1). The chlorination of PCB-169 is symmetric, and as such has no DM, but the chlorination pattern is such that all active sites (para and meta) are chlorinated and, thus can transfer electrons to, or take electrons from a neighboring receptor. In this case, toxic potency was determined by only reactivity of chlorinated sites. Because of its active surface, PCB-169 exhibits toxic potency that is lesser than that of congener 126, but greater than that of other congeners. Since these sites have already been established to be the most active, PCB-81 was classified as the third most potent congener. Finally PCB-77, which has a DM directed exactly perpendicular to the center axis does not preferentially expose para sites. However, it does expose the 3 and 3′ meta sites, which are both chlorinated. Due to this DM direction, PCB-77 is classified as being less potent than the 3 most potent congeners, but because it is non-ortho it is still more potent than the other 8 dioxin-like PCBs.

All mono-ortho PCBs were ranked similarly with TEFs of 0.00003 (Van den Berg et al., 2006) (Table 1). PCB-123 (2′,3,4,4′,5-pentachlorobiphenyl) was predicted to be less potent than PCB-189 (2,3,3′,4,4′,5,5′-heptachlorobiphenyl) because the DM does not expose any meta sites, whereas PCB-189 does. An angle θ ~ 90° is least toxic, because it does not expose para sites. The model derived was able to rank this group without ambiguity. Due to its mono-meta nature, PCB-74 was ranked as the 13th most potent congener. The model produces a negative DRF value for PCB-74 which indicates that the electron accepting nature of the hydrogenated para and meta sites outweighs the electron donating tendency of the chlorinated sites. In the DRF model, a value of 0 would be the threshold for dioxin-like toxicity, suggesting that PCB-74 has lesser toxic potency. Overall, the ranking of toxic potencies of PCB congeners suggested both DRF model and TEF were consistent. In addition, the DRF values of PCB congeners obtained from this study were well correlated with AhR binding affinity based on in vitro test (r² = 0.63, p < 0.01) (Safe et al., 1985; Mekenyan et al., 1996), indicating that the model was well-directed to predict the toxicity of PCB congeners. However, the DRF values were not correlated to the lipophilicity of PCB molecules such as log K_ow values (Svendsgaard et al., 1997). It is indicated that the log K_ow values of chemicals are not less explained the toxic potency in itself, but the values could be related to the bioaccessibility and bioavailability.

Bioaccumulation Characteristics

Bioaccumulation of PCBs is a complex and multifaceted process since metabolism to hydroxylated PCBs is species-specific and dependent to the types of enzymes present (Kaminsky et al., 1981; WHO, 2003). As aforementioned, the HOMO-LUMO gap is one of the main parameters relating to the electron-mediated bonding. That is, if there is any increase in the HOMO-LUMO energy gap of a PCB molecule when absorbed into an organism, this molecule is stabilized by tight binding to the AhR and becomes more difficult to be metabolized. Bioaccumulation of chemicals is determined by their rates of depuration. In fact, hydroxylation makes molecules more polar and thus more likely to be excreted directly or as conjugates. Thus, monitoring of changes in HOMO-LUMO gap before and after exposure of animals to PCBs would allow prediction of bioaccumulation of the untransformed molecules.

Synchrotron-radiation X-ray spectroscopy was employed for empirical determination of HOMO-LUMO gaps. Near-Edge X-ray Absorption Fine Structure (NEXAFS) and X-ray Emission Spectroscopy (XES) probed unoccupied and occupied molecular orbitals of PCBs, respectively and determined the HOMO-LUMO energy gaps by superposing the C 1 s NEXAFS and C Kα XES spectra. To determine center peak locations corresponding to the HOMO and LUMO states more precisely, 2nd derivatives of spectra were used (Bazylewski et al., 2011). Spectroscopic measurements were obtained for four PCBs: PCB-101 (2,2′,4,5,5′-pentachlorobiphenyl), PCB-105 (2,3,3′,4,4′-pentachlorobiphenyl), PCB-118 (2,3′,4,4′,5-pentachlorobiphenyl), and PCB-153 (2,2′4, 4′,5,5′-hexachlorobiphenyl). Superposition of non-resonant C Kα XES and C 1 s NEXAFS spectra of PCB-118 after microsome exposure (a) and their second derivatives (b), is given (Figure 5) as a representative example. This also demonstrates how the HOMO-LUMO energy gap was determined. The HOMO and LUMO levels were determined to be at 281.94 and 285.12 eV, respectively, which results in a gap of 3.18 eV. Measurements were repeated for four PCB congeners before and after exposure (Figure 6). It should be noted that the measured HOMO-LUMO gap of PCB-118 is smaller than that obtained by ab-initio calculations (4.82 eV). This is not unexpected due to the approximated estimation of the exchange-correlation functional used in the DFT calculation. Despite this, there is a clear trend that the HOMO-LUMO energy gap increases for all PCBs after exposure.

Figure 5. XES and XAS measurements. (A) Superposition of XES and XAS measurements of PCB-118 (exposed) (B) 2nd Derivative Method to determine HOMO-LUMO gap. HOMO-LUMO gap is determined by difference between π and π* energies, in this case 3.178 eV.

Figure 6. HOMO-LUMO gaps of 4 PCBs studied. The trend is for band-gap to increase (toxicity to decrease) after exposure. PCB-118 changes least, indicating most bio-accumulative, also indicated is an experimental band-gap threshold of toxicity, similar to Figure 1.

Change in the HOMO-LUMO gap is different for each PCB. PCB-153 has a gap change of 1.2 eV, whereas PCB-118 shows a relatively small change of 0.2 eV. This result suggests that PCBs become more recalcitrant to chemical reaction after exposure to microsomes and that PCB-153 becomes more difficult to be metabolized and thus more bioaccumulative than other PCBs. Based on these results and the information from Section First Principles Model of Toxicity, it can be concluded that a small HOMO-LUMO gap is a necessary condition to estimate the toxic potency while the bioaccumulative behavior is closely associated with the change in the HOMO-LUMO gap rather than the magnitude of the gap (molecule's toxicity).

Concluding Remarks: Environmental Implications

A model, based on first principles, was developed to predict the toxic potency of AhR-active compounds, including PCBs. The primary and necessary criterion for toxicity was determined to be “band-gap.” This separates more potent dioxin-like PCBs from less potent or non-dioxin-like PCB congeners. Dioxin-like PCBs were then examined more closely to determine “site-specific reactivity.” This reactivity was shown to be greater for para, followed by meta, and finally ortho sites. Mono-ortho-substituted PCBs had, on average, lesser total reactivity than did non-ortho PCBs. This reactivity model is consistent with established purely structural models, and in fact explains the structure models in a more physical, non-empirical way. Ranking of toxic potency of the 12 known dioxin-like PCBs was accomplished by inclusion of a directional reactivity factor which allowed accurate prediction of TEF for and ranking of all 12 dioxin-like PCBs. Another possible dioxin-like PCB, PCB-74, was identified and ranked 13th. The model did not provide direct prediction of established TEF values since these are consensus values determined from a meta-analysis of a number of different endpoints for a wide range of species. Even though this model only ranks the dioxin-like PCBs, which mediate toxicity through the AhR, it could be used to predict relative toxicity of other AhR-active compounds. The methodology can also be applied to other receptor-mediated effects as long as the LBD of the receptor can be defined. The model is most powerful when the receptor structure, particularly the DM and position of likely attachment is known. At a minimum the model is useful for identifying other halogenated biphenyls, like poly-brominated and mixed chloro-bromo analogs of the biphenyls, dioxins, furans and naphthalenes as well as alkylated polycyclic aromatic hydrocarbons that can bind to the AhR.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We gratefully acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC) and Canada Foundation for Innovation. Research described in this paper was performed at the Canadian Light Source, which is supported by NSERC, the National Research Council Canada, the Canadian Institutes of Health Research, the Province of Saskatchewan, Western Economic Diversification Canada, and the University of Saskatchewan. This work was also supported by the project entitled “Development of Technology for CO₂ Marine Geological Storage” and “Oil Spill Environmental Impact Assessment and Environmental Restoration” funded by the Korean Ministry of Land, Transport, and Maritime Affairs given to Prof. Jong Seong Khim. Prof. John P. Giesy was supported by the Canada Research Chair program, a Visiting Distinguished Professorship in the Department of Biology and Chemistry and State Key Laboratory in Marine Pollution, City University of Hong Kong. He was also supported by the program of 2012 “High Level Foreign Experts” (#GDW20123200120) funded by the State Administration of Foreign Experts Affairs, P.R. China to Nanjing University and the Einstein Professor Program of the Chinese Academy of Sciences.

References