|Home | About | Journals | Submit | Contact Us | Français|
DNA polymerases play a crucial role in the cell cycle due to their involvement in genome replication and repair. Understanding the reaction mechanism by which these polymerases carry out their function can provide insights into these processes. Recently, the crystal structures of human DNA polymerase λ (Polλ) have been reported both for pre- and post- catalytic complexes (García-Díaz et al., DNA Repair, 3, 1333, 2007). Here we employ the pre-catalytic complex as a starting structure for the determination of the catalytic mechanism of Polλ using ab initio quantum mechanical/molecular mechanical methods. The reaction path has been calculated using Mg2+ and Mn2+ as the catalytic metals. In both cases the reaction proceeds through a two step mechanism where the 3′-OH of the primer sugar ring is deprotonated by one of the conserved Asp residues (D490) in the active site before the incorporation of the nucleotide to the nascent DNA chain. A significant charge transfer is observed between both metals and some residues in the active site as the reaction proceeds. The optimized reactant and product structures agree with the reported crystal structures. In addition, the calculated reaction barriers for both metals are close to experimentally estimated barriers. Energy decomposition analysis to explain individual residue contributions suggests that several amino acids surrounding the active site are important for catalysis. Some of these residues, including R420, R488 and E529, have been implicated in catalysis by previous mutagenesis experiments on the homologous residues on Polβ. Furthermore, Polλ residues R420 and E529 found to be important from the energy decomposition analysis, are homologous to residues R183 and E295 in Polβ, both of which are linked to cancer. In addition, residues R386, E391, K422 and K472 appear to have an important role in catalysis and could be a potential target for mutagenesis experiments. There is partial conservation of these residues across the Pol X family of DNA polymerases.
DNA based organisms must replicate their genomes accurately and also remove and replace damaged nucleotides to maintain genome stability. These are critical processes since errors in replication and/or repair can result in mutations, some of which can lead to disease or even death [1, 2, 3]. DNA polymerases are the proteins in charge of the processes involved in DNA replication and repair. They are classified into several distinct families based on differences in the primary structure of their catalytic subunits .
All transactions associated with the processes involving DNA polymerases depend on the reaction catalyzed by these enzymes. This reaction is based on a nucleophylic attack of the O3′ of the primer-terminal nucleotide on the Pα of the incoming nucleotide to form an O-P bond, with subsequent release of pyrophosphate (PPi). The mechanism involved in this transfer has been hypothesized to occur with the help of two metal ions [5, 6]. Experimental and theoretical studies on the reaction mechanism of several polymerases support this hypothesis [7, 8, 9, 10, 11, 12, 13].
Theoretical studies on the reaction mechanism of DNA polymerase β (Polβ) have suggested that the reaction is largely associative. One study using QM calculations on model systems  has proposed a direct proton transfer from the O3′ to an O atom on the Pα. In another study, a full QM/MM treatment on the same system suggests the proton transfers from the primer O3′ to an aspartate residue in the active site . This latter mechanism is similar to the one proposed for T7 DNA polymerase [16, 17]. Conversely, theoretical calculations on the mechanism of DNA polymerase IV from Sulfolobus solfataricus have proposed the proton transfer to happen through an ordered water . This pathway has also been proposed previously as a likely mechanism for Polβ [9, 19]. In all cases the phosphate-breaking step has been proposed to be associative-like, although the general base may differ.
In addition to the chemistry involved in the catalytic step, the question arises as to the nature of the metal. Recent studies on DNA polymerase λ (Polλ) have suggested a slight preference for Mn2+ over Mg2+ as the activating metal based on lower activation energy for the former metal . However, for the homologous Polβ such preference is not observed .
Polλ is a family X polymerase that fills short gaps during DNA repair. It has been implicated in both non-homologous end-joining (NHEJ) of double stranded DNA breaks and in base excision repair of damaged bases [22, 23, 24, 25, 26]. Polλ shares 35% amino acid identity with Polβ, another family X polymerase [27, 28, 29, 30]. Experimentally, it has been shown that Polβ under-goes a large conformational change involving subdomain motion needed for catalysis . On the other hand this large subdomain motion is not observed in Polλ [31, 32]. However, this conformational change appears not to be rate limiting and it is possible that the rate limiting step is the same (or similar) for both Polλ and Polβ.
Recently, several X-ray structures of human Polλ have been obtained . These structures include a pre-catalytic ternary complex (pdb id 2PFO) with the primer terminus nucleotide poised for an in-line attack of the ribose O3′ to the Pα of the incoming non-hydrolyzable dUPNPP, with a Mn2+ in the catalytic metal site and a Mg2+ in the “dNTP-binding” position. In addition, a post-catalytic structure (pdb id 2PFQ) was obtained where a dCTP has been covalently added to the primer terminus with both metal sites occupied by Mn2+.
In the present contribution, we take the reported pre-catalytic X-ray structure (2PFO) as the starting point for quantum mechanical/molecular mechanical (QM/MM) calculations [34, 35, 36]. Two sets of calculations are performed with either Mg2+ or Mn2+ as the divalent metal ions in the active site to investigate whether there is a preference and/or different mechanism in the activation process. This is the first such theoretical study that compares the catalysis of a DNA polymerase with two different metals (Mg2+ and Mn2+) in the active site. Our calculations provide further support for the hypothesized two metal mechanism and show that both metals experience a significant charge transfer to some residues in the active site during the reaction. In addition, an energy decomposition analysis provides insight into the catalytic role of amino acid residues surrounding the active site, several of which are conserved among family X polymerases .
Molecular dynamics (MD) [38, 39, 40] and QM/MM [41, 42, 43, 44, 45, 46] calculations have been performed to investigate the reaction mechanism of human Polλ. MD simulations were carried out with the AMBER9 suite of programs . All QM/MM calculations were performed with modified versions of Gaussian03 and TINKER for the QM/MM reaction path calculations [48, 49]. The QM/MM optimizations were performed with an iterative method as described in  and  using the pseudobond model for all boundary atoms [34, 50]. Reaction paths were obtained using a combination of the quadratic string method (QSM) , reaction coordinate driving (RCD)  and quadratic synchronous transit (QST3) [53, 54].
The initial structure of Polλ was taken from the pre-catalytic X-ray structure (pdb id 2PFO) . This pre-catalytic structure contains a non-hydrolyzable nucleotide analog (dUPNPP), a Mg2+ as the dNTP-binding metal and a Mn2+ in the catalytic site. This structure was modified by replacing the dUPNPP by dUTP and the Mn2+ by Mg2+. Hydrogens were added using the XLEaP program in AMBER9 . This program was also used to solvate the structure in a box of 5856 water molecules. The system was minimized and subjected to molecular dynamics (MD) simulations using the PMEMD module of AMBER9 . A total of 2 ns of MD was carried out to equilibrate the system and the last snapshot was selected as the starting point for the QM/MM calculations. The final configuration was modified by retaining all atoms in the protein and all water molecules within 30 Å of the catalytic metal for the subsequent QM/MM calculations (see Supplementary Information for details).
Following equilibration, the QM subsystem was chosen to include both active site metals, side chains of D427, D429, D490, the primer dC nucleotide (excluding C5′ and the phosphate group), a part of the incoming dUTP nucleotide (triphosphate and C5′) and two water molecules that complete the coordination spheres of the metals, for a total of 72 QM atoms including 5 boundary atoms (Cα of aspartates, C5′ of dC and C4′ of dUTP). The total charge of the dUTP was set to −3 with one protonated oxygen on the γP as done in previous theoretical studies on Pol β and HIV reverse transcriptase [11, 15]. Note that we have not changed the incoming nucleotide to the canonical Watson-Crick pair (dTTP). This was done in order to avoid the docking of dTTP in the place of dUTP, and to remain as close as possible to the initial X-ray structure.
Two sets of calculations were performed with either Mg2+ or Mn2+ as active site metals. In the case of Mn2+ as the catalytic metal, all calculations were performed assuming a high spin state with all d electrons on both Mn atoms unpaired (see Supplementary Information) [55, 56, 57]. The QM subsystem was treated at the B3LYP level [58, 59] with a combined basis set where the 6–31G* basis was used for all reactive atoms (OD1 from D490, OD1 from D279, O in the ordered water, O3′, O3α and Pα), except for the Pα where an additional diffuse function was included. All other (non-reactive) QM atoms were represented with 3–21G. The Mn2+ atoms were treated with the LANL2DZ pseudopotential. The remaining atoms were included in the MM subsystem and treated with the parm99 force field .
Although there is a post-catalytic structure for Polλ , it contains a different nucleotide (dTTP) as the newly added base than the reactant structures. Because of this, the product structures for both systems (Mg2+ and Mn2+) were produced in silico from the reactant structure and subjected to QM/MM optimization. The optimized end-points were employed for the reaction path calculation with the quadratic string method (QSM) . This is a “chain-of-replica” method where the path is represented by a discrete number of structures [44, 60, 61, 62, 63]. In QSM the intermediate structures are obtained by a linear interpolation between the end points. Subsequently, all structures are optimized to the minimum energy path simultaneously. This affords an unbiased optimization of the path, which is highly efficient because all images on the path are calculated in parallel. The initial QSM paths consist of 10 points including end-points and were iteratively optimized with constrained MM-optimization [64, 65, 63]. After the initial paths were converged, the highest energy points were optimized to the closest transition state (TS) with the quadratic synchronous transit (QST3) method [53, 54]. Following the TS optimization the paths between the critical points were calculated with either QSM or RCD, see supplementary materials for details.
To verify the applicability of the basis set, an initial QSM path with 10 total points was calculated for the Mg2+ with all atoms being treated at the B3LYP/6-31G* level (6-31+G* for the Pα). This path was found to be the same as the ones calculated with the combined basis set and the barrier was ≈2.5 kcal/mol higher. In this case the highest energy point has not been optimized to the TS and the energy could drop around 1–1.5 kcal/mol. Therefore, based on this result we decided to perform the remaining calculations with the combined basis set for computational efficiency.
In this subsection we provide the results for the reaction mechanisms. In this case, we have decided to explore different reaction pathways based on previous studies [15, 18]. The different schemes representing the possible pathways differ in the general base used to deprotonate the O3′ of the primer nucleotide. These possible bases include D490, D429 or an ordered water that completes the coordination sphere of the catalytic metal (see Fig. 1).
The first scheme corresponds to the proton transfer from O3′ to the closest OD of D490. This scheme is equivalent to the one proposed by Lin et al. for human DNA polymerase β . In scheme 2, the proton is transferred to the closest OD in D429. Scheme 3 is based on a recently proposed mechanism for DNA polymerase IV from Sulfobolus solfataricus where the proton is transferred to an ordered water molecule . In this case only the structure with Mg2+ was taken into account to determine the proton acceptor.
The calculated potential energy barrier associated with the reaction path for scheme 1 is 17.6 kcal/mol (see Fig. S1 in supplementary materials ). The calculated reaction barrier for scheme 2 is very close to that of scheme 1, however, the optimized path (Fig. S2 in supplementary materials) shows that the reaction proceeds by an initial proton transfer to D490. This is followed by the breaking of the Pα-Oβ bond and ending with the direct proton transfer from D490 to D429. Based on this, we have ruled out scheme 2 since the initial step for this mechanism is equivalent to scheme 1.
In the case of the proton transfer to the ordered water, we obtained a reaction barrier of 35.4 kcal/mol. Another related mechanism involves a proton transfer through this water molecule to an O on the Pα . The reaction energy for this mechanism was found to be only 2 kcal/mol lower than the one for scheme 3. Note that in the present study, one oxygen is protonated in the Pγ, whereas Wang et al. used a deprotonated triphosphate for their calculations . In order to test this difference in protonation states, the reactant and product for scheme 3 were re-optimized with the deprotonated triphosphate. In our hands, the calculated reaction energy for scheme 3 in Polλ with the deprotonated triphosphate is 41.7 kcal/mol. In addition, scheme 1 was also tested with the deprotonated phosphate structures, with a resulting reaction energy of 25.2 kcal/mol. Experimentally, two proton transfers have been observed for this reaction , one corresponds to the deprotonation of the O3′ and the other is hypothesized to be the protonation of the phosphate. This protonation could be done by a protonated side chain or from a proton in the solvent. In this case, present and previous results suggest that the O3′ deprotonation is rate limiting, and could suggest that the protonation of the phosphate may happen after binding in the active site (as assumed in the present case), or after PPi formation as in [9, 14] with a lower barrier than for the first proton transfer. Also, experimentally it is not known if the PPi is deprotonated based on its pKa values which range between 7.0–9.0 .
Based on the above observations for the three different reaction mechanisms, our results point to scheme 1 with the protonated phosphate as the preferred mechanism for Polλ. This scheme is similar to the previously proposed mechanism by Lin et al.  for human Polβ where the proton is transfered to the homologous D256 residue. Conversely, in our hands the water-mediated mechanism proposed by Wang et al.  appears to not be viable for Polλ because of the high reaction energies associated with this scheme. From these results as well as those previously published it is difficult to generalize these reaction mechanisms to all the members of the X and Y family polymerases at this point [16, 17, 19, 14, 15, 9]. However, although difficult to generalize based on the general base, it appears that the mechanism is associative-like .
The relaxed and optimized reactant systems with Mg2+ and Mn2+ both show an RMSD of 1.7 Å for backbone atoms with respect to the 2PFO X-ray reference. This difference comes from the MD relaxation of the entire enzyme. However, the superposition of the active site residues (D427, D429, D490, dC, dUTP and metals) yields an RMSD below 0.3 Å for all heavy atoms in both Mg2+ and Mn2+ systems (see Figure 2). In the case of the product systems, the RMSD between the relaxed and optimized structures to the X-ray reference (2PFQ) is 1.9 Åfor all heavy atoms. Again, as in the reactant case, the RMSD for all heavy atoms in the active site residues remains around 0.2 Å (see Fig. 3). Note that this is a very good agreement, especially since the product structure was directly generated from the reactant via an optimization with the QM/MM method.
As seen in Figure 4 both the Mg2+ and Mn2+ paths show a large initial energy increase, which corresponds to the proton transfer. This is followed by a high energy, rugged valley characterized by the breaking of the Pα-Oβ bond. Finally the nucleotide is transferred and pyrophosphate is formed in the product structure. Note that in both cases, the product structures are around 4–5 kcal/mol higher in energy than the reactants. Overall, this mechanism is similar to the one proposed for Polβ, with the exception that the reaction energy for Polβ is negative and we do not observe a stable low energy intermediate after the proton transfer .
Two TSs are observed along the paths for Mg2+ and Mn2+ In the case of the Mg2+ the optimized transition states (TSs) give a total energy barrier of 17.6 and 16.4 kcal/mol with respect to the reactant and a frequency calculation revealed a single imaginary (negative) frequency for each maximum. For the Mn2+ the TSs are located at 16.8 and 17.6 kcal/mol respectively. A single imaginary frequency was obtained for each of these structures as well. Additionally in both cases the first TS was observed to correspond to the proton transfer and the second to the Pα-Oβ bond break (see Fig. 5).
The calculated barriers are close to the experimental counterpart, where a kcat of 6.0 s−1 has been reported for the correct insertion of dT opposite dA by human Polλ in the presence of Mg2+. This corresponds to a barrier of 16.6 kcal/mol , calculated from pre-steady state kinetic studies for the incorporation of a single nucleotide via the transition state formula (TST) . The barriers estimated using TST assume that the catalytic step is rate limiting and thus that the activation free energy ΔG‡ corresponds to the energy barrier. Note that our calculated barriers correspond to potential energy and the estimated experimental barriers correspond to free energy. Usually, the calculated potential energy barriers provide an upper bound to the corresponding free energy barriers in QM/MM calculations [36, 18]. Therefore, it can be expected that a free energy calculation would lower the value of the barriers closer to the experimental estimate.
Furthermore, it has been estimated experimentally that the energy barrier for the Mn2+ catalyzed reaction in human Polλ is between 0.4–0.9 kcal/mol lower than the corresponding Mg2+ one . In the present work, the energy difference in the barriers is not large enough to determine an ion preference or whether either TS is rate limiting. In addition, based on the approximations inherent in the methodology, the estimated experimental difference is too small and falls within the error of the method.
Figure 5 shows the superposition of the optimized transition states for both metals. As can be seen, the first TS of the Mg2+ system is very similar to the first TS of the Mn2+ system, where the proton has been transferred to D490. In the case of the second TS both structures show that the Pα-Oβ bond is broken. The only other major difference between the second TS of the Mg2+ and that of the Mn2+ systems is a slightly larger distance between the transfered proton and O3′ in the Mn2+ system. These trends are best observed in Figure 6, where selected distance changes along the reaction coordinate are presented. Note that as the reaction progresses from the first to the second TS, the O3′-Pα decreases considerably while the Pα-Oβ increases. Additionally, as can be seen in Figure 5, the Pα in the second TS shows a distorted trigonal bipyramidal phosphorane, resulting in an associative-like mechanism.
It is important to note that in our calculations the distance between the metals, either Mg2+ or Mn2+, remains almost constant along the path. The observed fluctuations for the inter-metal distances are less than 0.1 Å over the range of the reaction, with an average distance of 3.5 Å for Mg2+ and 3.7 Å for Mn2+, metal-ligand distances for reactants, TSs and products are given in the supplementary materials.
The analysis of the charge distribution in the QM subsystem presents an interesting picture of the reaction mechanism. Changes in charge for selected groups in the active site at the stationary points on the path (TSs and end-points) are shown in Table 1. As can be seen the overall charge transfer to D490, O3′, αP and the PPi are similar for the Mg2+ and Mn2+ catalyzed reactions, with a significant reduction in the negative charge on O3′ and a similar reduction in the positive charge on the PPi. Note however, that in the case of the Mg2+ both metals experience a large charge transfer to D427, D490 and the water molecule coordinated to the dNTP-binding Mg2+. On the other hand, for Mn2+ there is an observed charge transfer from the metals in the TSs, however, there is no overall change between the reactant and the product. The difference in charge transfer from the metals to the ligands may be due to the use of different basis sets/methodology for Mg2+ and Mn2+.
Further insights into the catalytic mechanism of Polλ may be obtained by understanding the role of individual residues in catalysis. To that end, the non-bonded interaction energy between the MM environment and the QM subsystem can be decomposed into contributions per residue [68, 69]. Note that this “breakdown” is only obtained in an approximate manner and is only qualitative. In this analysis the differences in interaction energies are calculated between individual residues and the QM subsystem when the system goes from reactant to TS(s). A negative contribution denotes a stabilizing contribution to the TS and a positive one denotes a destabilizing contribution. We consider a residue to provide a significant contribution to (de)stabilization if the interaction is greater than 1 kcal/mol, see supplementary materials for details.
Several residues are observed to have a significant effect on the TS barrier (see Fig. 7). Coulomb and Van der Waals analysis show that a significant part of the stabilization provided by the MM environment comes from the solvent (see Supplementary Material). Overall most of the interaction differences observed are electrostatic. The only Van der Waals (VdW) interactions that were observed were for the TSs of Mg2+ where an MM water has a contribution of approximately 3 kcal/mol, and for the second TS of the Mn2+ reaction, where R488 and a water molecule have an interaction of around 1 kcal/mol (see Supplementary Material).
From Fig. 7 it can be seen that there are a number of residues that have a significant contribution to the stabilization of the TSs. In particular R386, E391, R420, K422, K472, R488, E529, D574 and W575 are observed for both Mg2+ and Mn2+. None of these residues have been mutated experimentally on Polλ to our knowledge. However, the homologous residues to R420, R488 and E529 in Polβ have been experimentally studied. Note that the comparison between our calculations and the experimental mutagenesis results is only qualitative since the determined energy change is only a first order approximation to the TS (de)stabilization and does not include structural or dielectric screening contributions.
This residue shows an important stabilizing contribution to the TSs for both catalytic metals of around −1 kcal/mol in Mg2+ and −2 kcal/mol for Mn2+. Mutation of the equivalent Arg-183 to Ala in Polβ produces a reduction greater than 99% in enzyme activity . Moreover, this mutation has been observed in some esophageal cancers . Castro et al. have suggested the possibility of homologous residues in other polymerases to act as a general acid that protonates the pyrophosphate in their hypothesized two proton transfer mechanism . Incidentally, R420 happens to be the corresponding residue for this proposed mechanism. The protonation of the triphosphate by this residue is currently under investigation.
Experimental studies have shown that the mutation of the homologous R254 in rat Polβ to Ala decreases the kcat from 51 min−1 to 1 min−1 for an experimental ΔΔG of 2.33 kcal/mol . From Figure 8 it can be seen that R488 forms a salt bridge with D490. In our analysis, this residue provides a significant destabilizing contribution of around 6 kcal/mol. The protonation of D490 reduces its acid strength and thereby reduces the strength of the R488 interaction. The resulting destabilization of the TS by R488 can be attributed to the tight coupling between these two residues. This destabilization does not take into account any structural changes that result from the experimental mutation of this residue. The experimentally observed reduction in kcat could be due to the fact that the role of R488 may also be structural, by orienting D490 for a proper geometry of the coordination sphere of the catalytic metal and its proper positioning to deprotonate the O3′. Indeed, salt bridges between R488 and D490 are observed in the X-ray structures with Na, Mg and Mn (pdbid 2PFN 2PFO and 2PFQ respectively) with average distance of 2.8 Å. Therefore, by mutating the homologous R254 to A in Polβ, the salt bridge is removed and the structural integrity of the catalytic aspartate could be compromised. However, a minor alteration in the D256-R254 salt bridge distance is noticed (2.95 Å in the 1BPY X-ray structure).
Our results show that this residue also results in a destabilizing contribution of around 1 kcal/mol for both metals. This may be the result of the interaction with the dC base. Experimentally, this residue has been linked to gastric cancer in Pol β [1, 2]. Mutation of the homologous residue in Polβ to a lysine inhibits the base excision repair mechanism .
To our knowledge these are the only residues that have been studied experimentally and correspond to any of the residues suggested to have a significant contribution to catalysis in our residue analysis. The remaining residues could provide good targets for experimental mutagenesis studies. In particular K472 is located at the end of a loop region that is five residues longer in Polλ than in Polβ. The role of this loop and the residues involved in this loop are currently under experimental investigation. Although the residues obtained from this analysis are catalytically important, they have not been included in the QM subsystem due to computational limitations and due to the fact that it has been previously shown that residues that show significant non-bonded interactions with the QM subsystem can be adequately modeled with the MM force field .
To further understand the catalytic role of these residues, we have performed a sequence alignment of Polλ with the sequences of other X family polymerases as shown in Figure 8 bottom (see Fig. S11 in supplementary materials for full alignment), as in a previous study . This alignment, performed with T-Coffee , shows that the equivalent residues for R420 and R488 in Polλ are totally conserved among X family polymerases. K472 is highly conserved except in Polµ where there is an arginine instead. The remaining residues are conserved to some degree, either by identity in a couple of polymerases, or by residue function, e.g., charge of the residue. These results suggest that the conserved residues could have a role in the enzyme function and non-conserved residues could be involved in the specificity and/or selectivity of the different polymerases.
The reaction mechanism for human Polλ has been determined by means of QM/MM calculations for two catalytic metals. Three reaction schemes were tested to determine the general base that deprotonates the O3′. Our results point to D490 as the general base, the other two schemes (D429 or H2O) appear to not be viable in Polλ. In the case of D429 as the proton acceptor, the proton is transfered first to D490. For the H2O mediated scheme, the calculated reaction energies were too high. After determining the scheme, the reaction path was calculated for two systems, Mg2+ and Mn2+. In both cases, the reaction mechanism proceeds in two steps. The first step corresponds to the abstraction of the proton on O3′ by D490, followed by the breaking of the Pα-Oβ bond. This mechanism agrees with the previously proposed mechanism for human Polβ. The calculated energy barriers are around 17 kcal/mol for the Mg2+ and Mn2+ catalyzed reactions. These barriers are relatively close to the experimental determination of 16.6 kcal/mol for the incorporation of dT opposite dA with Mg2+. For both Mg2+ and Mn2+ a significant charge transfer is observed from the catalytic and dNTP-binding metals to some residues in the active site.
Energy decomposition analysis to explain individual residue contributions to catalysis was carried out. This analysis suggests that there are several residues that have a significant effect on the TS barrier. Sequence alignment studies show that some of these residues are strictly conserved among different enzymes of the Polymerase X family, pointing to a functional role for these amino acids. The other residues that are not conserved may be involved in specific enzyme requirements due to their unique roles. Among the calculated residues from the decomposition, three have been studied experimentally, R420, R488 and E529, on the corresponding residues on Polβ and two of these have been observed in certain cancers. The remaining residues have not been studied experimentally and could be an interesting target for mutagenesis studies.
Detailed methods for initial structure, QM/MM structure optimization and path optimization are given, as well as QSM paths for selected reaction schemes, metal-ligand distances, residue analysis for all residues and waters, and sequence alignment of Polλ with other Family X polymerases.
This research was supported by the intramural research program of the NIH and NIEHS. Computing time from the NIEHS and UNC computing facilities are gratefully acknowledged. GAC would like to thank Drs. S.K. Burger, Z. Lu and Prof. W. Yang at Duke University for the QSM code and Prof. Y. Zhang at NYU for enlightening discussions. We would like to thank the reviewers as well as Drs. T. Darden, B. Beard, G. Muller at NIEHS for critical reading of this manuscript.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Conflict of Interest statement
The authors declare that there are no conflicts of interest.