Abstract
The introduction of linear energy correlations, which explicitly relate adsorption energies of reaction intermediates and activation energies in heterogeneous catalysis, has proven to be a key component in the computational search for new and promising catalysts. A simple linear approach to estimate activation energies still requires a significant computational effort. To simplify this process and at the same time incorporate the need for enhanced complexity of reaction intermediates, we generalize a recently proposed approach that evaluates transition state energies based entirely on bondorder conservation arguments. We show that similar variation of the local electronic structure along the reaction coordinate introduces a set of general functions that accurately defines the transition state energy and are transferable to other reactions with similar bonding nature. With such an approach, more complex reaction intermediates can be targeted with an insignificant increase in computational effort and without loss of accuracy.
Introduction
In the last decade, developments in density functional theory (DFT) methods have made it possible to effectively perform computerbased screening studies in the design of heterogeneous catalysts^{1,2,3,4,5,6,7,8,9,10}. The combination of explicit DFTcalculated adsorption energies (i.e., descriptors) and simple modeling methods, which describe effectively and accurately the correlation between these descriptors and complex reaction and activation energies, is an instrumental tool in this approach. Models that are first order in energy, socalled linear scaling methods, have shown to provide sufficient accuracy when evaluating trends in catalytic activity^{11}. However, recently it has been shown that the linearity is closely connected to the intrinsic bondorder of the intermediate, thus suggesting that while accurate for the energy of stable intermediates on a potential energy surface (PES), it will be less accurate in describing the energy of metastable intermediates such as transition states^{12,13}. Instead we have introduced a new method, “rscaling”, which explicitly takes variations in bondorder along the reaction coordinate into account. This directly incorporates structural variations of transition states into the energy correlations and should potentially be more accurate than typical linear approaches, such as Brønsted–Evans–Polanyitype relations^{14,15,16,17,18,19}. In addition, it also substantially reduces the computational cost when establishing transition state (TS) energy correlations^{13,20,21}.
The rscaling scheme is built entirely on the bondorder conservation principle, which relies on the hypothesis that any transient structure (defined by the reaction coordinate r) between the initial and final state of an elementary reaction step along the minimum energy reaction pathway obey the following linear relation:
where E^{t} (r, E^{d}) and E^{d} are the chemisorption energies of a transient intermediate at a specific reaction coordinate r and that of a reaction relevant descriptor species, respectively. The constraints imposed by this hypothesis uniquely define both γ(r) and ξ(r), and the TS energies as a function of E^{d} can be obtained by maximizing the potential energy E^{t} (r, E^{d}) given by Eq. (1) along the reaction coordinate.
According to the bondorder conservation principle, the dissociation of a covalent bond on a surface is associated with a relocation of the electrons into the new bonds established between the surface and the fragments of the adsorbate such that the normalized bond order of the fragments is conserved^{22,23,24,25,26}. The electronic details of this process are characterized by the distinct type of bond that is broken and the associated functional forms of γ(r) and ξ(r) in Eq. (1), related to the specific orbitals in the bond and their interaction with the surface. This suggests that the cleavage of chemical bonds exhibiting the same local electronic environment can be grouped into classes. The members in each class expose a similar electron redistribution, and therefore, will share the same functional form of γ(r) and ξ(r).
In this work, we will justify this concept and demonstrate that the functions γ(r) and ξ(r), as defined in Eq. (1), are transferable between reactions within the same class. This allows us to obtain the TS energyscaling relations for a large number of reactions from a database of template functions of γ(r) and ξ(r) determined accurately for a single element in the class. The functional form for all elements within the class can be easily generated by rescaling the template functions, which requires fewer calculations than establishing the functions for each individual reaction. In the following, we use the dehydrogenation reactions of sp^{3} hybridized C, N, and O atoms to apply this concept. We show explicitly that when using the functions γ(r) and ξ(r) derived from CH_{4}, NH_{3}*, and H_{2}O* bond activation as our template functions to define each class, we can effectively and accurately derive the TS energyscaling relations of any other sp^{3}hybridized C–H, N–H, and O–H dissociation reaction. This not only broadens our understanding of chemical reactions and how the reaction energy evolves, it also significantly speeds up our computational approach.
Results
Application of rscaling to dehydrogenation reactions
We have studied dehydrogenation reactions of hydrocarbons, ammonia, amine, water, and alcohols, which all contain sp^{3}hybridized bonds, to illustrate the underlying similarity between specific bondbreaking reactions and how this can be used to address a larger set of reactions at little extra computational cost. The reaction steps examined are listed below. By considering the activation of CH_{4}, NH_{3}*, and H_{2}O* molecules on fcc(211) surfaces as representatives for dehydrogenation reactions of species containing C, N, and O centers, we applied the rscaling scheme to derive the functional form of γ(r) and ξ(r) for each reaction and identify transient intermediate energies as a function of r and E^{d} through Eq. (1), as shown in Fig. 1.
For all transient structures with fixed bond lengths ranging from 1.2 to 2.2 Å, the slope and intercept of the scaling of chemisorption energies (i.e., E^{t} as a function of E^{d}) changes uniformly as the bond is splitting (Fig. 1a–c). The observed increase in the slope stems from the bondorder conservation and is a manifestation of the increasing surfaceadsorbate bond order as the R–H bond is breaking^{20,21}. The surface corresponding to E^{d} = 0 in Fig. 1g–i identifies the intercept ξ, which can be interpreted as the potential energy of the reaction proceeding on a noninteracting surface. On the basis of the linear relations in Eq. (1), the monotonically increasing functions γ(r) and ξ(r) can be fitted as shown in Fig. 1d–f and the 2D PES for the three reactions shown in Fig. 1g–i can be derived from Eq. (1). The first order saddle points on the 2D PES (blue lines) define the TS energies, which are obtained by maximizing Eq. (1) along the reaction coordinate r:
Equation (17) defines the mapping of these saddlepoint energies into a set of functions depending only on a single descriptor (E_{ads}(A))_{A=C,N,O}. These functions represent the TS energyscaling relations (solid black lines in Fig. 1g–i). When compared with the actual calculated TS energies, we find a mean absolute error (MAE) of 0.07, 0.09, and 0.08 eV for the activation of CH_{4}, NH_{3}*, and H_{2}O*, respectively, as shown in Fig. 1j–l. Using the rscaling scheme, TS scaling relations of all the dehydrogenation reactions can be derived as shown in Supplementary Figure 7.
The γ(r) function in Eq. (1) reflects the intrinsic variation in the bondorder between the surface and the reactant and therefore, it changes independently of the orientation of the catalyst surface. In contrast, the function ξ(r) in Eq. (1) has a very strong dependence on the surface orientation, and therefore needs to be determined specifically for a certain reaction coordinate. However, when bondbreaking processes involve a specific bonding network along a welldefined reaction coordinate, such as sp^{3}hybridized bonds (e.g., CH_{4}, *CH_{3}, *CH_{2}, and *CH), these processes will form a class of reactions where all members share similar functional forms of γ(r) and ξ(r). As a result, each reaction in the class will experience similar electron redistribution along the reaction coordinate as reflected in the comparison of molecular orbital and projected densities of states for C_{1} dehydrogenation reactions shown in Fig. 2a and b. In the following, we justify that through a simple rescaling, the functions γ(r) and ξ(r) can be transferred between bond dissociation reactions belonging to the same class. The corresponding PES can then be found through Eq. (1) and the TS energies can be derived from Eq. (17).
Once the members of a specific class have been identified, a successful rescaling of the functions γ(r) and ξ(r) within the class depends on only two conditions: (1) that the template functions γ_{ i }(r) and ξ_{ i }(r) derived explicitly for a single element i in the class are accurately established and (2) that the initial and final values of the functions for other members j ≠ i are determined. Rescaling the template functions according to the initial and final points of the jth member, allows us to generate the complete set of functions for all j, as illustrated in Fig. 3. In the following, we show how the transfer scheme can be applied to a class of reactions involving C–H bond splitting. The activation of CH_{4} is used to establish the set of template functions γ(r) and ξ(r) from which new functions are generated to describe the activation of CH_{3}*, CH_{2}*, and CH*. The domain for any set of the continuous functions γ(r) and ξ(r) is reduced to a generalized reaction coordinate involving only the bond distance, which in this work ranges from 1.2 to 2.2 Å. The functional values at r = 1.2 and 2.2 Å will be obtained from the linear relations described by Eq. (1) using only two metal surfaces as shown in Fig. 4a and b. The two metals are chosen, based on the established scaling relations for the equilibrium initial and final state structures of the elementary reaction step, such that the smallest energyaveraged rootmeansquare error (σ_{avg}) is achieved (see Methods)^{21}. For CH_{3}*, CH_{2}*, and CH* activation, the two metals with the smallest σ_{avg} are {Ag,Ru}, {Ag,Ru}, and {Au,Rh}, respectively. Figure 4a and b shows the relations for CH_{3}*, CH_{2}*, and CH* as established using the two metal surfaces for C–H bond lengths fixed at r = 1.2 and 2.2 Å, respectively. The obtained values for slopes γ and intercepts ξ are seen to be very close to those obtained from linear relations fitted using all nine metal surfaces (Table 1) and these values are used as rescaling parameters to perform the transfer process as shown in Fig. 4c and d. TS energies for CH_{3}*, CH_{2}*, and CH* activation can now be obtained using Eq. (17). The whole process is detailed in the Methods section.
Discussion
The general scheme and the electronic structure similarities of the bonds formed between the catalyst surface and the reactants allow us to expand upon the three reaction classes defined by the functions obtained from CH_{4}, NH_{3}*, and H_{2}O* activation. In principle, all reactants R–H, where R is any residue forming an sp^{3}hybridized bond with hydrogen through C, N, or O should fall in one of the three classes defined above, e.g., activation of CH_{3}*, CH_{2}*, CH*, CH_{3}CH_{2}*, and CH_{3}CH_{2}CH_{2}* belong to the same class as CH_{4} activation, NH_{2}*, NH*, and CH_{3}NH* to the same class as NH_{3}* activation, and OH*, CH_{3}OH*, C_{2}H_{5}OH*, and C_{3}H_{7}OH* to the same class as H_{2}O* activation. Using the transfer scheme, we established the γ(r) and ξ(r) functions for all the members of each class and subsequently extrapolated the TS energyscaling relations of all the considered dehydrogenation reactions from Eq. (17) (Supplementary Figure 8). In Fig. 5, a comparison is shown between the transition state energies obtained using this simple model and the energies obtained using a constrained optimization in DFT based on either the fixed bondlength method or the Nudged–Elastic–Band (NEB) approach^{27,28}. Our comparison shows an overall MAE of only 0.12 eV using CH_{4}, NH_{3}*, and H_{2}O* activation on (211) surfaces to obtain the template functions for each class. Using alternative members to establish the template functions, the overall MAE ranges from 0.11 to 0.17 eV. It is noteworthy to mention that our model provides a much cheaper route for obtaining TS energies, since only four constrained calculations are required to secure the necessary information for a specific reaction belonging to a class with preestablished template γ(r) and ξ(r) functions. In contrast to traditional NEB methods, where TS energies are computed for every reaction step and catalyst surface, our model classifies reactions and applies a transfer scheme to effortlessly derive the TS energies thus significantly reducing the associated computational cost.
In summary, the functions γ(r) and ξ(r) are shown to depend only on the electronic nature of the activated bonds irrespective of the surface type or residual groups. Using three series of dehydrogenation reactions involving sphybridized C, N, and O on three low index metal surfaces: (100), (111), and (211), we have shown that the derived functions γ(r) and ξ(r) define three unique classes of reactions and that they can be transferred easily between bondbreaking reactions belonging to the same class. Only four constrained calculations are required to apply the transfer scheme for a given reaction step, allowing access to the full PES from which TS energies along the reaction coordinate can be derived. Moving forward, this enables us to build databases of bond specific γ(r) and ξ(r) functions for a welldefined set of surface reaction classes and then employ the transfer scheme to obtain relevant kinetic parameters. This approach introduces a significant reduction in the computational cost of calculating transition state energies. A simple estimation shows a factor of ~36 in the reduction of cpu time compared to traditional NEB methods, resulting in significant savings in power consumption as well as time spent by scientific personnel (Methods). This expansion of scaling methods to include the description of reaction transients increases the efficiency in the computation of kinetic parameters without loss of accuracy and primarily it enables an even faster screening of catalyst materials.
Methods
DFT calculations
DFT calculations are performed using the planewavebased PWSCF (QuantumESPRESSO) code and the Atomic Simulation Environment (ASE)^{29}. The ultrasoft Vanderbilt pseudopotential method with BEEF exchangecorrelation functional is adopted^{30,31,32,33}. A cutoff energy of 500 eV for the wavefunctions is used. The fcc(211), fcc(111), and fcc(100) surfaces of nine transition metals including Cu, Ni, Ag, Pd, Re, Rh, Ru, Pt, and Au are chosen as catalysts for the dehydrogenation reactions. The (211) slabs are built with 4 atomic layers in the [111] direction in rectangular 1 × 3 supercells with the bottom three [111] layers fixed during structural optimizations. The (111) and (100) slabs are built with four atomic layers with the bottom two layers fixed during structural optimizations using hexagonal 3 × 3 and rectangular 2 × 2 supercells, respectively. The Monkhorst–Pack scheme is used for sampling the Brillouin zone^{34} and kpoint grids of 4 × 2 × 1, 4 × 4 × 1 and 4 × 4 × 1 are selected for (211), (111), and (100) slabs, respectively. The vacuum thickness is set to more than 10 Å in all slab calculations and the dipole correction is applied to decouple the interaction between slabs. During structural relaxations, the residual force between atoms is converged to a value below 0.02 eV/Å.
TS energies for the dehydrogenation reactions are calculated using the fixed bondlength (FBL) method as implemented in ASE. All TS structures show a single imaginary frequency thus identifying it as a first order saddlepoint. The fixed bondlength calculations range from 1.2 to 2.2 Å with increments of 0.1 Å, which is accurate enough to cover transition states of any dehydrogenation reaction. TS’s with a bondlength shorter than 1.2 Å are close to the initial state, and in these cases the reaction is treated as if it has no barrier. TS energy calculations for all reactions (2–16) are performed on fcc(211) surfaces and for reactions (3), (6), and (7) calculations have also been performed on fcc(111) and fcc(100) surfaces.
Fitting procedure
In the rscaling scheme, fitting of the slope γ(r) and the intercept ξ(r) is done applying the scipy.optmize.minimize method using a fifthorder polynomial. The functions are constrained to be monotonically increasing.
Choosing the two metal surfaces
Established linear scaling relations for the initial and final state equilibrium structures for a given reaction have been used for selecting the two transitionmetal surfaces. Selection of a metal surface is partly determined by the rootmean square error (σ) of the linear regressions of the initial and final state:
where R_{IS} and R_{FS} are the metalspecific absolute residuals of the initial and final state linear regression, respectively. The energy span of the descriptor is also a significant factor. To ensure that the metals chosen reflect both the reactive as well as the noble surfaces, we define an energyaveraged rootmean square error of two metals (m_{1} and m_{2}) as:
The choice of two transitionmetal surfaces that minimizes Eq. (21) is finally used for building the scaling relations of transient structures with bondlength fixed at r = 1.2 and 2.2 Å. More details about the method selecting the two metals can be found in ref. ^{21}.
Steps of transfer scheme for obtaining TS energy relations

(1)
Obtain the γ(r) and ξ(r) functions for a single reaction in a class using a global or simplified rscaling scheme^{21}. These functions will be used as the template functions for the reaction class.

(2)
For a given reaction in the class, calculate adsorption energies (E_{rel}) of the initial state and final state on a sufficiently large number of metal surfaces and build the linear relations of E_{rel} vs. descriptor for the two states, from which two metal surfaces are determined as detailed above. Note that the initial and final state structures in this step are the equilibrium structures, not those with bondlength fixed at 1.2 and 2.2 Å.

(3)
Calculate the adsorption energies (E_{rel}) on the two metal surfaces determined in step (2) of two transient structures with bondlength fixed at 1.2 and 2.2 Å, respectively (four calculations total disregarding the clean surfaces and gasphase calculations). Then build the linear relations of E_{rel} vs. descriptor for the two transient structures using just two points for each (Fig. 4a, b) and obtain the values of γ(1.2), ξ(1.2), γ(2.2), and ξ(2.2) from the linear relations.

(4)
Rescale the template functions γ(r) and ξ(r) functions as shown in Fig. 3.

(5)
Build the TS energyscaling relation for the reaction according to Eq. (17).
NEB versus rscaling comparison
To establish the transition state scaling relation for a specific reaction using nine catalyst surfaces, the NEB method requires:
The rscaling method needs:
where N_{core} is the number of cpucores used per image or structure. Considering that NEB calculations take longer time to converge, we shall assume that NEB calculations take a factor of two more time than the constrained minimization calculations used in the rscaling method. The difference in the number of cores used is (2×72N_{core})/4N_{core} = 36, in favor of rscaling. Thus rscaling introduces a substantial reduction in cpu time, which renders it a highly efficient and accurate method for estimating TS energies.
Data availability
Part of the data sets generated and analysed during the current study are included in the Supplementary Information file and others are available in the “figshare” repository, https://doi.org/10.6084/m9.figshare.5601622.v2.
References
 1.
Honkala, K. et al. Ammonia synthesis from firstprinciples calculations. Science 307, 555–558 (2005).
 2.
Van Santen, R. A. & Neurock, M. Concepts in theoretical heterogeneous catalytic reactivity. Catal. Rev. Sci. Eng. 37, 557–698 (1995).
 3.
Norskov, J. K., Bligaard, T., Rossmeisl, J. & Christensen, C. H. Towards the computational design of solid catalysts. Nat. Chem. 1, 37–46 (2009).
 4.
Greeley, J., Jaramillo, T. F., Bonde, J., Chorkendorff, I. & Norskov, J. K. Computational highthroughput screening of electrocatalytic materials for hydrogen evolution. Nat. Mater. 5, 909–913 (2006).
 5.
Greeley, J. & Mavrikakis, M. Alloy catalysts designed from first principles. Nat. Mater. 3, 810–815 (2004).
 6.
Greeley, J. et al. Alloys of platinum and early transition metals as oxygen reduction electrocatalysts. Nat. Chem. 1, 552–556 (2009).
 7.
Linic, S. & Christopher, P. Overcoming limitation in the design of selective solid catalysts by manipulating shape and size of catalytic particles: epoxidation reactions on silver. ChemCatChem. 2, 1061–1063 (2010).
 8.
Salciccioli, M. & Vlachos, D. G. Kinetic modeling of Pt catalyzed and computationdriven catalyst discovery for ethylene glycol decomposition. ACS Catal. 1, 1246–1256 (2011).
 9.
Mpourmpakis, G. & Vlachos, D. G. Computationalbased catalyst design for thermochemical transformations. MRS Bull. 36, 211–215 (2011).
 10.
Montemore, M. M. & Medlin, J. W. Scaling relations between adsorption energies for computational screening and design of catalysts. Catal. Sci. Technol. 4, 3748–3761 (2014).
 11.
AbildPedersen, F. et al. Scaling properties of adsorption energies for hydrogencontaining molecules on transitionmetal surfaces. Phys. Rev. Lett. 99, 016105 (2007).
 12.
CalleVallejo, F., Loffreda, D., KoperMarc, T. M. & Sautet, P. Introducing structural sensitivity into adsorption–energy scaling relations by means of coordination numbers. Nat. Chem. 7, 403–410 (2015).
 13.
AbildPedersen, F. Computational catalyst screening: scaling, bondorder and catalysis. Catal. Today 272, 6–13 (2016).
 14.
Bell, R. P. The theory of reactions involving proton transfers. Proc. R. Soc. Lond. A 154, 414–429 (1936).
 15.
Evans, M. G. & Polanyi, M. Inertia and driving force of chemical reactions. Trans. Faraday Soc. 34, 11–24 (1938).
 16.
Cheng, J. et al. Brønsted−Evans−Polanyi relation of multistep reactions and volcano curve in heterogeneous catalysis. J. Phys. Chem. C 112, 1308–1311 (2008).
 17.
Bligaard, T. et al The Brønsted–Evans–Polanyi relation and the volcano curve in heterogeneous catalysis. J. Catal. 224, 206–217 (2004).
 18.
Wang, S. et al. Universal Brønsted–Evans–Polanyi relations for C–C, C–O, C–N, N–O, N–N, and O–O dissociation reactions. Catal. Lett. 141, 370–373 (2011).
 19.
van Santen, R. A. & Neurock, M. Theory of surface chemical reactivity. Russ. J. Phys. Chem. B 1, 261–291 (2007).
 20.
Plessow, P. N. & AbildPedersen, F. Examining the linearity of transition state scaling relations. J. Phys. Chem. C. 119, 10448–10453 (2015).
 21.
Yu, L. & AbildPedersen, F. Bond order conservation strategies in catalysis applied to the NH3 decomposition reaction. ACS Catal. 7, 864–871 (2017).
 22.
Shustorovich, E. Chemisorption phenomena: analytic modeling based on perturbation theory and bondorder conservation. Surf. Sci. Rep. 6, 1–63 (1986).
 23.
Shustorovich, E. Heat of molecular chemisorption from bondorderconservation viewpoint: why morse potentials are so efficient. Surf. Sci. 181, L205–L213 (1987).
 24.
Shustorovich, E. Metal effects in the Fischer–Tropsch synthesis: bondorderconservationmorsepotential approach. Catal. Lett. 7, 107–118 (1990).
 25.
Shustorovich, E. & Sellers, H. The UBIQEP method: a practical theoretical approach to understanding chemistry on transition metal surfaces. Surf. Sci. Rep. 31, 1–119 (1998).
 26.
van Santen, R. A. On Shustorovich’s bondorder conservation method as applied to chemisorption. Recl. Trav. Chim. PaysBas 109, 59–63 (1990).
 27.
Henkelman, G. & Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 113, 9978–9985 (2000).
 28.
Henkelman, G., Uberuaga, B. P. & Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 113, 9901–9904 (2000).
 29.
Paolo, G. et al. QUANTUM ESPRESSO: a modular and opensource software project for quantum simulations of materials. J. Phys. Condens. Matter 21, 395502 (2009).
 30.
Vanderbilt, D. Soft selfconsistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B 41, 7892–7895 (1990).
 31.
Laasonen, K., Pasquarello, A., Car, R., Lee, C. & Vanderbilt, D. Car–Parrinello molecular dynamics with Vanderbilt ultrasoft pseudopotentials. Phys. Rev. B 47, 10142–10153 (1993).
 32.
Laasonen, K., Car, R., Lee, C. & Vanderbilt, D. Implementation of ultrasoft pseudopotentials in ab initio molecular dynamics. Phys. Rev. B 43, 6796–6799 (1991).
 33.
Wellendorff, J. et al. Density functionals for surface science: exchangecorrelation model development with Bayesian error estimation. Phys. Rev. B 85, 235149 (2012).
 34.
Monkhorst, H. J. & Pack, J. D. Special points for Brillouinzone integrations. Phys. Rev. B 13, 5188–5192 (1976).
Acknowledgements
We acknowledge support from the U.S. Department of Energy, Office of Basic Energy Sciences, to the SUNCAT Center for Interface Science and Catalysis at SLAC National Accelerator Laboratory and Stanford University.
Author information
Affiliations
Contributions
L.Y. conducted all DFT calculations for the N–H and O–H bonds dissociation reactions and processed and analyzed data. L.V. conducted all DFT calculations for C–H bond dissociation reactions. F.A.P. conceived the calculations and analyzed results. All the authors interpreted the structures and contributed to the preparation of the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing financial interests.
Additional information
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Yu, L., Vilella, L. & AbildPedersen, F. Generic approach to access barriers in dehydrogenation reactions. Commun Chem 1, 2 (2018). https://doi.org/10.1038/s420040170001z
Received:
Accepted:
Published:
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.