Numerical calculation of protein-ligand binding rates through solution of the Smoluchowski equation using smoothed particle hydrodynamics
- Wenxiao Pan^{1},
- Michael Daily^{2} and
- Nathan A Baker^{3}Email author
DOI: 10.1186/s13628-015-0021-y
© Pan et al.; licensee BioMed Central. 2015
Received: 17 January 2015
Accepted: 30 March 2015
Published: 7 May 2015
Abstract
Background
The calculation of diffusion-controlled ligand binding rates is important for understanding enzyme mechanisms as well as designing enzyme inhibitors.
Methods
We demonstrate the accuracy and effectiveness of a Lagrangian particle-based method, smoothed particle hydrodynamics (SPH), to study diffusion in biomolecular systems by numerically solving the time-dependent Smoluchowski equation for continuum diffusion. Unlike previous studies, a reactive Robin boundary condition (BC), rather than the absolute absorbing (Dirichlet) BC, is considered on the reactive boundaries. This new BC treatment allows for the analysis of enzymes with “imperfect” reaction rates.
Results
The numerical method is first verified in simple systems and then applied to the calculation of ligand binding to a mouse acetylcholinesterase (mAChE) monomer. Rates for inhibitor binding to mAChE are calculated at various ionic strengths and compared with experiment and other numerical methods. We find that imposition of the Robin BC improves agreement between calculated and experimental reaction rates.
Conclusions
Although this initial application focuses on a single monomer system, our new method provides a framework to explore broader applications of SPH in larger-scale biomolecular complexes by taking advantage of its Lagrangian particle-based nature.
Keywords
Diffusion Smoluchowski equation Smoothed particle hydrodynamics Protein-ligand interactions Binding rates AcetylcholinesteraseBackground
In the “perfect” enzyme [1] acetylcholinesterase (AChE), the rate-limiting step for catalysis is diffusional encounter [2,3]. Specifically, the active site lies at the bottom of a 20 Å-deep gorge, and the diffusion of substrate into it is accelerated by electrostatic steering [4,5]. Its diffusion-limited behavior, complex geometry, and strong electrostatic influence has made AChE a useful target for both experimental and computational studies of biomolecular diffusion [4-11].
Two major classes of methods have been used to estimate diffusion rates in biomolecular systems. Mesoscopic coarse-grained methods like Monte Carlo [12-14], Brownian dynamics (BD) [8,9,15], and Langevin dynamics [16,17] simulations trace the trajectories of individual coarse-grained particles driven by Brownian motion. Such simulations typically consider dilute ligand concentrations so that electrostatic protein-ligand interactions can be modeled by the Poisson-Boltzmann equation [18,19] with a few notable exceptions [20]. Alternatively, continuum models can be used to treat the diffusion of ligand concentration in space around a biomolecule by the Smoluchowski equation [6,21-25]. In particular, an adaptive finite element approach [26] has been used to numerically solve the Smoluchowski equation, and it shows higher accuracy in predicting experimental data about the ligand binding rates than the coarse-grained BD modeling [6]. For dilute ligand concentrations, electrostatic interactions can also be modeled with the Poisson-Boltzmann equation like the mesoscale approach [6,7]. However, for more concentrated ligand solutions, continuum models can also model the electrostatic potential near the biomolecular surface using a regularized Poisson-Nernst-Planck formulation [24,25], allowing screening of the ligand-receptor interactions by its time-dependent distribution around the protein.
Here, we follow the continuum approach but solve the Smoluchowski equation using a new smoothed particle hydrodynamics (SPH) method [27,28]. Unlike Eulerian grid-based methods such as finite element method (FEM), SPH is a Lagrangian particle-based method. SPH has been used with good accuracy for numerically solving partial differential equations (PDEs) describing momentum, mass and energy conservation laws [27]. In SPH, the domain is discretized into a set of “particles” that serve as interpolation points to numerically solve the governing PDEs. The SPH discretization of PDEs is based on a meshless interpolation scheme, which allows the PDEs to be written in the form of a system of ordinary differential equations (ODEs). SPH has a straightforward discretization without the need for time-consuming FEM mesh construction around complicated geometries such as biomolecules. Due to its Lagrangian nature, SPH has many advantages for modeling physical phenomena involving moving boundaries, large deformation of materials, multiphases, and advection-dominated diffusive transport [28-30]. Specifically, in SPH, free surfaces and interfaces between fluids move with particles, and hence, there is no need for front tracking schemes. And the non-linear advection term is embedded in the material derivative in the Lagrangian coordinate system, and hence, SPH models advection exactly. In addition, the similarity of SPH to molecular dynamics and mesoscopic coarse-grained particle methods (e.g., dissipative particle dynamics, BD, and Langevin dynamics), allows coupling of simulations across scales to build a multiscale modeling framework. This is our primary goal with the current work: to enable the multiscale and multiphysics description of biomolecular dynamics and ligand recognition. To the best of our knowledge, SPH has not been widely used in modeling biomolecular systems. Thus, in the present work, we aim to take the first step to introduce SPH into this field through the development of a SPH model for biomolecular diffusion with AChE as a test case.
In the SPH model, the Smoluchowski equation is numerically solved and the ligand binding rates are calculated from flux across the reactive boundary as in the previous studies using FEM [6,21-25]. However, in the previous FEM studies, active sites were modeled using the absolute absorbing (Dirichlet) boundary condition (BC). This BC has a simple description on the reactive boundaries but assumes infinitely fast chemical reactions between the enzyme and the ligand; i.e., a “perfect enzyme”. In our model, we take into account imperfect and non-instantaneous reactivity and thus solve the equation subject to a reactive (Robin) BC.
To solve the Smoluchowski equation subject to Robin BC using SPH, we use a continuum surface reaction method [31] which we have recently adapted to solve the Navier-Stokes equations subject to slip (Robin) boundary conditions [32]. In this formulation, the Robin BC is replaced by a reflective Neumann BC and a source term added into the governing equation. The derivation of the method is based on the approximation of the sharp boundary with a diffuse interface of finite thickness by means of a color function. This method is general for any arbitrary complex geometries and thus appropriate for modeling Robin BC in biomolecular systems with complex structures.
Results and discussion
Spherical test systems
Before the numerical method was applied to a biomolecular system with complicated geometry, we verified it on simple spherical test cases. Specifically, we considered a diffusing sphere with a radius R _{1}. The entire domain was confined by the outer boundary Γ _{ b } determined as a spherical surface with the radius of R _{2}=125 Å.
Application to AChE-ligand binding rates
We applied the SPH method to study the ligand binding kinetics of a simple spherical cationic ligand to the mouse acetylcholinesterase (mAChE) under various ionic strength conditions. Specifically, we performed the time-dependent calculations at ionic strengths of 0.0, 0.05, 0.10, 0.15, 0.20, 0.50 and 0.67 M until the diffusion reaches the steady-state. To achieve the highest accuracy with affordable computational cost, a resolution of Δ x=2 Å was used in all the following calculations.
Comparison of Debye-Hückel fits vs. ionic strength between experiment and simulations
Data | \({k_{\text {on}}^{0}}\) | \({k_{\text {on}}^{\mathrm {H}}}\) | Z _{ E } | RMSD _{ exp } |
---|---|---|---|---|
Experimental [4] | 9.8±0.6 | 1.30 | 2.3±0.2 | 0.00 |
BD | 9.1±0.3 | 2.66 | 1.7±0.2 | 1.52 |
SMOL FEM [7] | 9.5±0.1 | 1.92 | 2.8±0.1 | 0.37 |
SPH (Dirichlet) | 10.2±0.1 | 2.19 | 2.9±0.1 | 0.57 |
SPH (Robin, α=8×10^{3}) | 9.6±0.1 | 1.88 | 3.0±0.1 | 0.33 |
Hexa-mutant (experiment) | 1.8±0.1 | 0.57 | 1.2±0.2 | 0.00 |
Hexa-mutant (SPH, Robin) | 2.23±0.00 | 1.77 | 2.37±0.02 | 0.88 |
We also assessed the accuracy of SPH method for describing the ligand-binding kinetics of a mAChE surface mutant. We tested the surface hexa-mutant (E84Q, E91Q, D280V, D283N, E292Q, and D372N) from Radic et al. [4], which reduces the reaction rate by about a factor of 4 across the 0 to 0.67 M ionic strengths. For the mutants, which are nearly isosteric with the wild-type protein, we used the same SPH model as the wild type, but recalculated the electrostatic potentials for the mutant charge distribution. As presented in Table 1, the Robin BC SPH model has qualitative accuracy: predicting \(k_{\text {on}}^{0}\) of 2.23 compared to 1.80 from Radic et al. [4] with a k _{on} RMSD of 0.88 over the entire ionic strength range studied.
Conclusions
The Robin BC offers a new way to incorporate reactive surfaces into continuum diffusion models for rate calculations. This Robin-based model incorporates a new parameter α, which has units of Å/ μs and can be related to the probability of reaction within distance Δ x to the boundary and time interval Δ t by \(P=1-\exp (-\alpha \frac {\Delta t}{\Delta x})\) [33]. Thus, α=0 corresponds to zero reactivity (reflective Neumann BC) while α=∞ corresponds to absolute reactivity (absorbing Dirichlet BC).
There are two possible origins for the differences between the current SPH model results and past FEM calculations using the Dirichlet BC. First, the current SPH work uses a more recent mAChE structure (4B82) while the previous FEM calculations used an older structure (1MAH). Second, our SPH model uses a fixed resolution uniformly on both solution domain and boundaries, while the FEM adaptively meshes the reactive boundary with higher resolution.
This work has provided an initial demonstration that the Lagrangian (particle-based) SPH method out-performs the Eulerian (grid-based) FEM [6] in accurately predicting ligand binding rates in AChE. This result is important because while both methods can be used to study molecules of the size of AChE, SPH is more scalable to larger systems such as the synapse geometry where AChE operates. Additionally, due to its Lagrangian nature, SPH can easily incorporate other physical phenomena such as fluid flow or protein flexibility.
We have demonstrated that superior performance can be achieved using a probabilistic reactive (Robin) BC rather than a simple Dirichlet BC. In fact, the Robin BC is likely more biologically relevant than the Dirichlet BC. While the AChE enzyme is considered nearly “perfect” with a diffusion-limited reaction rate, there is experimental evidence that a very small fraction of substrates entering the active site gorge do not react. Specifically, recent kinetic experiments suggest that through unknown mechanisms, the PAS limits the rate of progression of non-substrates of any size to the catalytic site [34]. In addition, molecular dynamics simulations suggest that the PAS provides a selective gating function, for example by fluctuations in the gorge width that are likely to let acetylcholine but not let larger molecule pass through [35,36].
Methods
Governing equation and boundary conditions
specifying a Dirichlet BC on the outer boundary Γ _{ b } where the concentration is equal to a bulk concentration p _{bulk}. The outer boundary is often a spherical surface with a radius chosen to ensure that the ligand-protein potential is spherically symmetric and/or can be approximated analytically [6]. For the current study with mAChE, this outer boundary has radius R _{2}≈128 as determined following a procedure similar to Song et al. and Chen et al. [6,23]. Also following Song et al. and Chen et al., p is normalized such that p _{bulk}=1.
SPH discretization of equations and boundary conditions
In this section, we present SPH discretization of the Smoluchowski equation, using Eq. 2 if the Dirichlet BC is used and Eq. 9 if the Robin BC is assumed. To simplify notation, we omit superscript r for the variables in Eq. 9 in the subsequent derivations.
where q=|r|/h. With this form of weighting function, only particles within 2h distance from particle i contribute to the summations in the SPH equations. We have chosen h=1.3Δ x where Δ x is the size of the cubic lattice.
In the simulations presented below, we set h _{ r }=h but it could be set differently in future applications.
Calculation of potentials of mean force
We calculated the potential of mean force W(x) using the recently published 2.1 Å resolution structure of mAChE [39]. To prepare this structure for the calculation, we assigned titration states of ionizable residues using PROPKA [40] at pH 7, and we used PDB2PQR v1.8 [41,42] to assign atomic radii and charges. APBS v1.4 was used to perform a nonlinear Poisson-Boltzmann multi-grid calculation of the electrostatic potential over the entire simulation domain [43]. The small and large domains were set to 600 Å and 400 Å, respectively, with a fine grid spacing of 0.600 Å. For APBS calculations, we used the single Debye-Hückel boundary condition, a smoothed molecular surface, and protein and solvent dielectrics of 2 and 78.54, respectively. Atomic charges were mapped onto the grids using cubic B-spline discretization. The calculated potential was mapped onto the SPH discretization points of protein and solvent via trilinear interpolation.
Appendix A: continuum surface reaction method
subject to the reflective Neumann BC (Eq. 10). Note that ϕ is non-zero only in Ω _{ a }, where it is equal to 1 as defined in Eq. 28. Thus, the modified governing equation takes its final form as Eq. 9.
Appendix B: calculation of reaction rate
Declarations
Acknowledgments
The authors gratefully acknowledge the funding support by the Applied Mathematics Program within the U.S. Department of Energy Office of Advanced Scientific Computing Research as part of the Collaboration on Mathematics for Mesoscopic Modeling of Materials (CM4). The research was performed using PNNL Institutional Computing at Pacific Northwest National Laboratory. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract DE-AC06-76RL01830. This research was also supported by National Institutes of Health grant GM069702 to NAB.
Authors’ Affiliations
References
- Eigen, M. 1974. Diffusion control in biochemical reactions. In Quantum statistical mechanics in the natural sciencespp, ed. Mintz SL and Widemeyer SM. 37–61. New York: Plenum.View ArticleGoogle Scholar
- Bazelyansky, M, E Robey, and JF Kirsch. 1986. Fractional diffusion-limited component of reactions catalyzed by acetylcholinesterase. Biochem.25(1): 125–30. doi:10.1021/bi00349a019.View ArticleGoogle Scholar
- Nolte, HJ, TL Rosenberry, and E Neumann. 1980. Effective charge on acetylcholinesterase active-sites determined from the ionic-strength dependence of association rate constants with cationic ligands. Biochem.19(16): 3705–11. doi:10.1021/bi00557a011.View ArticleGoogle Scholar
- Radic, Z, PD Kirchhoff, DM Quinn, JA McCammon, and P Taylor. 1997. Electrostatic influence on the kinetics of ligand binding to acetylcholinesterase - Distinctions between active center ligands and fasciculin. J Biol Chem.272(37): 23265–77. doi:10.1074/jbc.272.37.23265.View ArticleGoogle Scholar
- Tan, RC, TN Truong, JA McCammon, and JL Sussman. 1993. Acetylcholinesterase - electrostatic steering increases the rate of ligand-binding. Biochem.32(2): 401–3. doi:10.1021/bi00053a003.View ArticleGoogle Scholar
- Song, YH, YJ Zhang, TY Shen, CL Bajaj, A McCammon, and NA Baker. 2004. Finite element solution of the steady-state Smoluchowski equation for rate constant calculations. Biophys J.86(4): 2017–29.View ArticleADSGoogle Scholar
- Song, YH, YJ Zhang, CL Bajaj, and NA Baker. 2004. Continuum diffusion reaction rate calculations of wild-type and mutant mouse acetylcholinesterase: Adaptive finite element analysis. Biophys J.87(3): 1558–66. doi:10.1529/biophysj.104.041517.View ArticleADSGoogle Scholar
- Wade, RC, ME Davis, BA Luty, JD Madura, and JA McCammon. 1993. Gating of the active-site of triose phosphate isomerase - Brownian dynamics simulations of flexible peptide loops in the enzyme. Biophys J.64(1): 9–15.View ArticleADSGoogle Scholar
- Tara, S, AH Elcock, PD Kirchhoff, JM Briggs, Z Radic, P Taylor, and al et. 1998. Rapid binding of a cationic active site inhibitor to wild type and mutant mouse acetylcholinesterase: Brownian dynamics simulation including diffusion in the active site gorge. Biopolymers.46(7): 465–74.View ArticleGoogle Scholar
- Davis, ME, JD Madura, BA Luty, and JA McCammon. 1991. Electrostatics and diffusion of molecules in solution - Simulations with the University-of-Houston-Brownian dynamics program. Comput Phys Commun.62(2–3): 187–97.View ArticleADSGoogle Scholar
- Davis, ME, JD Madura, J Sines, BA Luty, SA Allison, and JA McCammon. 1991. Diffusion-controlled enzymatic-reactions. Methods Enzymology.202: 473–97.View ArticleGoogle Scholar
- Genest, D. 1989. A Monte-Carlo simulation study of the influence of internal motions on the molecular-conformation deduced from two-dimensional NMR experiments. Biopolym.28(11): 1903–11. doi:10.1002/bip.360281107.View ArticleGoogle Scholar
- Saxton, MJ. 1992. Lateral diffusion and aggregation - a Monte-Carlo study. Biophys J.61(1): 119–28.View ArticleADSGoogle Scholar
- Kerr, RA, TM Bartol, B Kaminsky, M Dittrich, J-CJ Chang, SB Baden, and al et. 2008. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM J Sci Comput.30(6): 3126–49.View ArticleMATHMathSciNetGoogle Scholar
- Northrup, SH, JO Boles, and JCL Reynolds. 1988. Brownian dynamics of cytochrome-c and cytochrome-c peroxidase association. Sci.241(4861): 67–70. doi:10.1126/science.2838904.View ArticleADSGoogle Scholar
- Eastman, P, and S Doniach. 1998. Multiple time step diffusive Langevin dynamics for proteins. Proteins-Structure Funct Genet.30(3): 215–27.View ArticleGoogle Scholar
- Yeomans-Reyna, L, and M Medina-Noyola. 2001. Self-consistent generalized Langevin equation for colloid dynamics. Phys Rev E. 64(6). Article No. 066114. doi:06611410.1103/PhysRevE.64.066114.
- Baker, NA.2005. Biomolecular applications of Poisson-Boltzmann methods. Reviews in Computational Chemistry(KB Lipkowitz, R Larter, and TR Cundari, eds.), Vol. 21. New York: Wiley-VCH, Inc.Google Scholar
- Baker, NA. 2005. Improving implicit solvent simulations: a Poisson-centric view. Curr Opinion Struct Biol.15(2): 137–43.View ArticleGoogle Scholar
- Senapati, S, CF Wong, and JA McCammon. 2004. Finite concentration effects on diffusion-controlled reactions. J Chem Phys.121(16): 7896–900.View ArticleADSGoogle Scholar
- Kurnikova, MG, RD Coalson, P Graf, and A Nitzan. 1999. A lattice relaxation algorithm for three-dimensional Poisson-Nernst-Planck theory with application to ion transport through the gramicidin A channel. Biophys J.76(2): 642–56. doi:10.1016/s0006-3495(99)77232-2.View ArticleGoogle Scholar
- Smart, JL, and JA McCammon. 1998. Analysis of synaptic transmission in the neuromuscular junction using a continuum finite element model. Biophys J.75(4): 1679–88.View ArticleGoogle Scholar
- Cheng, Y, JK Suen, D Zhang, SD Bond, Y Zhang, Y Song, and et al.2007. Finite element analysis of the time-dependent Smoluchowski equation for acetylcholinesterase reaction rate calculations. Biophys J.92(10): 3397–406. doi:10.1529/biophysj.106.102533.View ArticleADSGoogle Scholar
- Lu, B, MJ Holst, JA McCammon, and YC Zhou. 2010. Poisson-Nernst-Planck equations for simulating biomolecular diffusion-reaction processes I: Finite element solutions. J Comput Phys.229(19): 6979–94. doi:10.1016/j.jcp.2010.05.035.View ArticleADSMATHMathSciNetGoogle Scholar
- Lu, B, and JA McCammon. 2010. Kinetics of diffusion-controlled enzymatic reactions with charged substrates. PMC Biophys.3: 1. doi:10.1186/1757-5036-3-1.View ArticleGoogle Scholar
- Holst, MJ. 2001. Adaptive numerical treatment of elliptic systems on manifolds. Adv Comput Mathematics.15(1): 139–91. doi:10.1023/a:1014246117321.View ArticleMATHMathSciNetGoogle Scholar
- Monaghan, JJ. 2005. Smoothed particle hydrodynamics. Rep Prog Phys.68: 1703–59.View ArticleADSMathSciNetGoogle Scholar
- Monaghan, JJ. 2012. Smoothed particle hydrodynamics and its diverse applications. Ann Rev Fluid Mech.44: 323–46.View ArticleADSMathSciNetGoogle Scholar
- Pan, W, AM Tartakovsky, and JJ Monaghan. 2013. Smoothed particle hydrodynamics non-Newtonian model for ice sheet and ice shelf dynamics. J Comput Phys.242(1): 828–42.View ArticleADSMATHMathSciNetGoogle Scholar
- Pan, W, D Li, AM Tartakovsky, S Ahzi, M Khraisheh, and M Khaleel. 2013. A new smoothed particle hydrodynamics non-Newtonian model for friction stir welding: Process modeling and simulation of microstructure evolution in a magnesium alloy. Int J Plasticity.48: 189–204.View ArticleGoogle Scholar
- Ryan, EM, AM Tartakovsky, and C Amon. 2010. A novel method for modeling Neumann and Robin boundary conditions in smoothed particle hydrodynamics. Comput Phys Commun.181: 2008–23.View ArticleADSMATHMathSciNetGoogle Scholar
- Pan, W, J Bao, and AM Tartakovsky. 2014. Smoothed particle hydrodynamics continuous boundary force method for Navier-Stokes equations subject to a Robin boundary condition. J Comput Phys.259: 242–59.View ArticleADSMathSciNetGoogle Scholar
- Gillespie, DT. 2007. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem.58: 35–55.View ArticleADSGoogle Scholar
- Beri, V, JT Auletta, GM Maharvi, JF Wood, AH Fauq, and TL Rosenberry. 2013. Hydrolysis of low concentrations of the acetylthiocholine analogs acetyl(homo)thiocholine and acetyl(nor)thiocholine by acetylcholinesterase may be limited by selective gating at the enzyme peripheral site. Chem Biol Interactions.203(1): 38–43. doi:10.1016/j.cbi.2012.09.017.View ArticleGoogle Scholar
- Zhou, HX, ST Wlodek, and JA McCammon. 1998. Conformation gating as a mechanism for enzyme specificity. Proc Natl Acad Sci U S A.95(16): 9280–3.View ArticleADSGoogle Scholar
- Baker, NA, and JA McCammon. 1999. Non-Boltzmann rate distributions in stochastically gated reactions. J Phys Chem B.103(4): 615–7.View ArticleGoogle Scholar
- Zhou, HX. 1990. On the calculation of diffusive reaction rates using Brownian dynamics simulations. J Chem Phys.92(5): 3092–95.View ArticleADSGoogle Scholar
- Allen, M, and D Tildesley.1989. Computer Simulation of Liquids. Oxford: Clarendon Press.Google Scholar
- Andersson, CD, N Forsgren, C Akfur, A Allgardsson, L Berg, C Engdahl, and et al.2013. Divergent structure-activity relationships of structurally similar acetylcholinesterase inhibitors. J Med Chem.56(19): 7615–24. doi:10.1021/jm400990p.View ArticleGoogle Scholar
- Li, H, AD Robertson, and JH Jensen. 2005. Very fast empirical prediction and rationalization of protein p K _{ a } values. Proteins Structure Funct Bioinformatics.61(4): 704–21. doi:10.1002/prot.20660.View ArticleGoogle Scholar
- Dolinsky, TJ, P Czodrowski, H Li, JE Nielsen, JH Jensen, G Klebe, and et al.2007. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res.35: 522–5. doi:10.1093/nar/gkm276.View ArticleGoogle Scholar
- Dolinsky, TJ, JE Nielsen, JA McCammon, and NA Baker. 2004. PDB2PQR: an automated pipeline for the setup of poisson-boltzmann electrostatics calculations. Nucleic Acids Res.32: 665–7. doi:10.1093/nar/gkh381.View ArticleGoogle Scholar
- Baker, NA, D Sept, S Joseph, MJ Holst, and JA McCammon. 2001. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc Natl Acad Sci U S A.98(18): 10037–41. doi:10.1073/pnas.181342398.View ArticleADSGoogle Scholar
- Lange, R-J. 2012. Potential theory, path integrals and the Laplacian of the indicator. J High Energy Phys.2012(11): 1–46.View ArticleADSGoogle Scholar
Copyright
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.