Exploration of nonlinear parallel heterogeneous reaction pathways through Bayesian variable selection

Inversion is a key method for extracting nonlinear dynamics governed by heterogeneous reaction that occur in parallel in the natural sciences. Therefore, in this study, we propose a Bayesian statistical framework to determine the active reaction pathways using only the noisy observable spatial distribution of the solid phase. In this method, active reaction pathways were explored using a Widely Applicable Bayesian Information Criterion (WBIC), which is used to select models within the framework of Bayesian inference. Plausible reaction mechanisms were determined by maximizing the posterior distribution. This conditional probability is obtained through Markov chain Monte Carlo simulations. The efficiency of the proposed method is then determined using simulated spatial data of the solid phase. The results show that active reaction pathways can be identified from the redundant candidates of reaction pathways. After these redundant reaction pathways were excluded, the controlling factor of the reaction dynamics was estimated with high accuracy.


Introduction
Nonlinear reaction dynamics in natural sciences are governed by reaction mechanism variables, such as reaction rate constants, diffusion coefficients, and reaction pathways [1,2]. Therefore, understanding the nonlinear dynamics of reaction and mechanism are important for scientific theory and ensuring precise simulations in Earth sciences. Many of the important reactions in natural sciences are heterogeneous, i.e., the reaction occurs at the interface between two or more phases, which are typically minerals and water [1,3]. Understanding the mechanism of heterogeneous reactions is important for multiple disciplines, including environmental science, geochemistry, hydrogeology, oceanography, nuclear waste disposal, and soil science. Inversion is a common way of understanding the mechanism of heterogeneous reactions; however, the complexity and nonlinearity of heterogeneous reactions, which typically involve sequences of several reactions occurring in parallel, makes inversion difficult.
Many laboratory experiments have attempted to understand the mechanism of heterogeneous reactions. These experiments typically involve quantitative observations of either (1) the temporal evolution of the concentrations of aqueous species [4][5][6][7], or (2) the spatial distribution of minerals or concentrations of aqueous a e-mail: royanagi@jamstec.go.jp (corresponding author) species at a single time [8][9][10][11]. The reaction rate constant in a reaction pathway can be determined inversely from the former by analyzing the temporal evolution of the solution chemistry and equilibration. However, it is impossible to determine the reaction rate and reaction pathways from the latter because the stability of a mineral cannot be determined without fluid composition data. Therefore, it is important to develop a versatile method that can extract nonlinear dynamics directly from a spatial dataset. The number of model parameters to be estimated is an important factor for inversion analysis. In heterogeneous parallel reactions, many possible reaction pathways are present, and many reaction rate constants need to be estimated. However, in general, the extent of fitting increases with the number of model parameters used for fitting. This phenomenon is well known as overfitting, which is a very common concept in the field of machine learning [12][13][14]. When overfitting occurs, the estimated model parameters are substantially affected by noise, which can reduce the predictive ability of the model. Thus, to estimate reliable reaction rate constants by inversion, reaction rate constants of unimportant reaction pathways must be identified and excluded from redundant candidates of reaction pathways.
In this study, we propose a framework for exploring the reaction pathways of nonlinear and parallel heterogeneous reactions from the spatial distribution of minerals alone. We used the concept of Bayesian estima-tion from machine learning [12,13]. The concept has been applied to various fields including physics [15][16][17][18][19][20], brain science [21], astronomy [22], and Earth sciences [23]. Notably, the Bayesian estimation has been applied to identify reaction pathways in the fields of biology and chemical engineering [24][25][26][27][28][29][30]. In these fields, the estimation of reaction rates from noisy observations is the technically important and challenging topic. Biochemical reaction networks are characterized by nonlinear systems that can be changed by spatio-temporal scale. These characteristics are similar to reactions in the Earth sciences where heterogeneous reactions between rock and solution are coupled with diffusive and/or advective transport of solution. Therefore, the Bayesian estimation can be an effective way to identify heterogeneous reaction pathways in the Earth sciences.
Under Bayesian estimation, a statistical model can be selected by evaluating the Bayesian free energy, and the important variables can be identified among redundant variables from a set of candidate models and given data [12]. Typically, huge computational costs are required to estimate the Bayesian free energy. Therefore, the Schwarz information criterion [31] that approximates Bayesian free energy has been used, but its usage is validated only if the statistical model is regular (i.e., the posterior distribution can be approximated by normal distribution). Recently, an alternative approach that gives an approximate value of the Bayesian free energy and can be used for both regular and singular statistical models has been proposed [32]. This approach has made it significantly easier to apply Bayes' theorem for nonlinear model selection problems. Using the framework proposed in this study, we successfully extract important reaction pathways of heterogeneous reactions using only the observable spatial distribution of minerals after the reaction. This paper is organized as follows. In Sect. 2, we introduce the method for estimating the reaction rate constants form the observed spatial mineral distribution. We first present a forward modeling approach for an observed spatial distribution of solids. Bayesian inference is employed to construct the posterior distribution, i.e., the conditional probability of the model parameters for an observed spatial distribution of phases. In Sect. 3, we validate the proposed estimation method by determining the reaction mechanism from a synthetic spatial distribution of solid phases. Section 4 presents the discussion and summary.

Method for estimating Model parameters
In this study, we analyze the spatio-temporal dynamics of a nonlinear heterogeneous reaction. To reflect prospected experimental constraints, we suppose a situation in which temporal changes are not observable and spatial changes at a single time are observable. Figure 1 shows a conceptual model of parallel heterogeneous reactions. In this conceptual model, three minerals (M a , M b , and M c ) are present at the initial time in discretized time t and space x for a finite difference approach. At each grid point, a set of partial differential equations conserving the masses of the chemical components in the system is solved, considering diffusion and heterogeneous reactions.y Heterogeneous mineral-water reactions are included in the model by using reaction rates that depend on the solution composition at each x and t (Cx,t). Active reaction pathways are determined from the local variable Cx,t. At the observable time t = T , the type and amount of mineral present varies with x because the reaction pathways differ according to Cx,t t = 0. These minerals can react and either be consumed or produced by each other. For example, one reaction can be written in the following generalized form: where A (aq) is an aqueous species in the solution, and k ab and k ba (s −1 ) are the rate coefficients of the reaction pathway for M b after M a and M a after M b , respectively. Other reactions can be written similarly; in this heterogeneous reaction network, the total number of reaction pathways L is given by L = N m (N m − 1) where N m is the total number of minerals. Notably, the value of k l is non-zero when the reaction pathway l is active but zero when the reaction pathway l is inactive. This concept is used for identifying important reaction pathways, as described later. Whether the reaction pathways are active depends on the local variable of concentration of A (aq) (C (mol cm −3 ); Fig. 1). The concentration C varies in space x and time t according to the reaction and the diffusion of A (aq) . The amount of mineral i (m (i) (mol cm −3 )), varies with x and t because it is a function of C.
To reflect realistic experimental environments, C is set as an unobservable value. The observable values in this study are the spatial variation of the amount of minerals and the volume fraction of the solution (porosity) at the start (t = 0 in Fig. 1) and end of the reaction (t = T in Fig. 1). Note that discretized space and time step (x, t) are used in Fig. 1 rather than continuous space and time (x , t )

Forward modeling of mineral-fluid interactions using a reaction-diffusion model
The spatio-temporal evolution of C x ,t can be modeled using the mass-conservation equation of aqueous species [1]; where ϕ (−) is porosity, D (cm 2 /s) is the constant of diffusivity for aqueous species, and Ξ l (C) (mol cm −3 ) is the loss or gain density of aqueous species according to the reaction pathway l. The first term on the right side of Eq. (1) expresses the concentration change due to diffusion of the aqueous species, whereas the second term expresses the change due to heterogeneous reactions. The rate of the heterogeneous reaction in a reaction pathway l, ∂Ξ l (C)/∂t (mol cm −3 s −1 ), can be written as a linear equation as a function of concentration, as follows: where C l (eq) (mol cm −3 ) and k l (s −1 ) are the constant of equilibrium concentration and the rate constant for a reaction pathway l, respectively. From the rate of heterogeneous reactions, the rate of increase or decrease in the amount of mineral i can be written as where m (i) (mol cm −3 ) and v l i are the amount of mineral i and stoichiometric coefficient of mineral i relative to aqueous species in the reaction pathway l, respectively. ϕ can be calculated from one minus the sum of the volume fraction of all minerals, as follows: where V i (cm 3 mol −1 ) is a constant representing the molar volume of mineral i. When Eqs. (1)-(3) are solved, x and t need to be discretized. The discretization of Eq. (1) is described in the Appendix. As a result of discretization, the concentration (C x,t ), amount of mineral i (m (i) x,t ), and porosity (ϕ x,t ) for a discretized space x (∈ {1, . . . , X}) and time t (∈ {1, . . . , T }) are obtained through forward modeling by inputting the following unknown parameter vector Θ : For a given parameter set Θ , the theoretical amount of all observation series (i.e., minerals and porosity) at the observable target time T with dimensions of J = N m +1 is The output for observation series j, y , is the sum of the response function for the corresponding input y (j) x and the observation noise x ; that is, At this point, we assume that the observation noise follows a Gaussian distribution with a mean of zero and a standard deviation of σ j : for x ∈ X obs , where X obs is a set of observable spatial points. Thus, in the forward modeling, the conditional probability of the observed spatial distribution of the mineral and porosity using the model parameters Θ is written as } is the set of observed series at position x. This formula can be used to obtain the probabilistic estimation of the mineral distribution after reaction as the conditional probability for the given model parameters. This can then be compared with experimental data.

Exhaustive search method for exploring reaction pathways
This section describes the method for achieving the first goal of this study, i.e., to identify the active reaction pathways. After the active reaction pathways are identified, the value of the rate constants are estimated (goal 2) in Section II C.
We can determine whether a reaction pathway l is active or inactive using the k l value of a reaction pathway l. Thus, Θ can be formulated as where k l is the intrinsic rate constant of the reaction pathways, the symbol • represents the Hadamard product, and Θ = {k 1 , . . . , k L }. c is an L dimensional binary vector, defined as Each variable c l is 0 or 1; c l = 1 if the lth variable belongs to the combination and is used for inversion. In contrast, c l = 0 if it does not . The term c is used to represent the indicators (Fig. 2). This formulation is key for clarifying one aspect of the program setting; i.e., identifying the active reaction pathways. The best c for modeling an objective variable is determined by minimizing the Bayesian free energy F using the exhaustive search (ES) method, which is the simplest variable selection method [33,34]. In the ES method, all combinations of used and unused variables are exhaustively searched, which requires estimation of 2 L combinations of variable. With the ES method, the F(c) is calculated for all combinations of explanatory variables; the combination that minimizes the F(c) is determined as the optimal combination.
Using Bayes' theorem, the posterior probability can be expressed as proportional to the likelihood function and the prior probability, as follows: Here, Y X = {y ex 1 , . . . , y ex X } and X is the number of spatial data points in the single observation series. We assume a continuous uniform prior probability p(c), where the posterior distribution is proportional to a marginalized likelihood function defined as follows: In addition, we assume that previous quantitative measurements of the minerals and porosity do not affect the present measurement; thus, Because reaction rate constants must be positive and may vary by orders of magnitude, as for p(Θ|c), we assume that the case of c l = 1, p(k l |c l = 1) has a continuous uniform distribution in logarithmic space, as follows: are the maximum and minimum values of the rate coefficient for the reaction pathway l within the range of searched values. In the case of c l = 0, p(k l |c l = 0) is assumed to be the delta function, as follows: The negative logarithm of the marginalized likelihood is termed the Bayesian free energy F, which is defined as The minimization of F is identical to the posterior probability maximization. However, it is difficult to analytically calculate the integral over the parameter Θ in the marginal likelihood. A well-established solution is thermodynamic integration [16] that numerically calculates F using the Markov chain Monte Carlo (MCMC) method [35]; however, this method generates huge computational costs because the expectations over several distributions with different pseudo-temperatures must be calculated. In this study, we numerically calculate the WBIC for the Bayesian model selection of reaction pathways. The WBIC value gives an approximate value of F [32]. This approach is expected to substantially reduce the computational costs of numerical calculation. The WBIC is defined as where n is the total number of data points (n = XJ) and β is the inverse pseudo-temperature (β > 0). L n is a negative logarithmic likelihood function, defined as As defined in Eq. (18), WBIC is identical to the average nL n over the posterior distribution with β = 1/ ln(n), whose β is different from the standard Bayesian estimation of the posterior (β = 1). The MCMC method [36] is employed to obtain WBIC values. In the MCMC simulations, the candidate is generated using a Gaussian proposal density (N (0, 0.01)). The same proposal is used for all models. The candidate is accepted with the Metropolis-type probability from the transition from Θ to Θ * as follows: where E is an energy function, Using these equations, WBIC values for 2 7 combinations of c were calculated using the MCMC method with the setting of β = 1/ln(n). The vector of active reaction pathways, i.e., the reaction pathways with nonzero rate constants, can be given by the indicator vector that minimizes the WBIC value.

Estimation of model parameters
Afterĉ is identified, the parameter in Θ can be estimated using Bayes' theorem, as follows: The likelihood function p(Y X |Θ,ĉ) is then given by Thus, the posterior distribution can be expressed as Fig. 3 Schematic illustration of the reactive transport model considered in this study [37]. At x = 0, SiO 2(aq) is transported by diffusion at a rate of D through the porous media of powdered olivine, and secondary minerals (talc, serpentine, and brucite) are formed. The gray bold line shows the overall reaction between the two minerals. k l is the effective rate constant of the reaction pathway l which is obtained by the MCMC method.

An example of parallel heterogeneous reactions
In this section, we validate the proposed method using simulated data. As an example of the potential applications of the proposed methodology, we consider coupled diffusion and reaction in the presence of olivine, quartz, and H 2 O (Fig. 3). This is chosen for the following reasons: (1) The diffusional metasomatic zoning of talc and serpentine zones between quartz and olivine has been classically modeled [38][39][40]; (2) this heterogeneous reaction network is relatively simple, but involves the dissolution/precipitation processes of several minerals as well as element diffusion [8,41].

Validation dataset
The black dotted line in Fig. 4 (a-e) shows the mineral spatial distribution, which is obtained by Θ = {10 −3.301 , 10 −3.301 , 0, 0, 10 −3.301 , 0, 0}; i.e., k 1 , k 2 , and k 5 are non-zero. The equilibrium concentrations are set as {C 1 (eq) ,. . . , C 7 (eq) } = {10 −6.57 , 10 −9.41 , 10 −5.86 , 10 −5.86 , 10 −7.54 , 10 −8.01 , 10 −8.01 } and D is set to 10 −4.2 [41]. The temperature-pressure variable V i at 300 • C and 10 MPa was obtained from the petrological software Perple_X [42], which can determine physical properties of minerals thermodynami-cally. We assume that the mineral spatial distribution can be observed at steps of 0.03 cm and is represented by adding Gaussian noise. To impose the same level of Gaussian noise on each observable series, σ j was assumed to be σ j = σ × y After adding Gaussian noise, we set negative values of the synthesized spatial data to zero, because the amount of mineral cannot be negative. The circle in Fig. 4(ae) shows the synthesized mineral spatial distribution. Each observed series contains 60 data points (X = 60) and there are five observation series (J = 5); thus, the total number of data points is 300 (n = 300). The parameters for prior information k (min) l and k (max) l are set to 10 −8 and 10 0 , respectively. The noisy dataset is then analyzed using the proposed method to determine its effectiveness for extracting and estimating the model parameters Θ.

Validation result
10 4 samples were obtained using the MCMC method with the setting of β = 1/ ln(n), 1000 burn-in time, and 10 thinning intervals, which were used to calculate WBIC values for 2 L candidates of c. As shown in a later section (Sect. 3.5), autocorrelation of chains is generally low in each MCMC run, suggesting that chains are well mixed. Figure 5 shows the calculated value of WBIC for  (Fig. 5a-d). In contrast, another 32 models have the highest WBIC values, and these worst models are uniformly ranked at 97 (Fig.  5a-d). The ranking result shows that c 1 , c 2 , and c 5 = 1 appeared frequently in combinations ranked at < 15 ( Fig. 5a, b). In contrast, c 3 , c 4 , c 6 , and c 7 = 0 frequently ranked at < 15 ( Fig. 5a, b), indicating that k 3 , k 4 , k 6 , and k 7 are comparably unimportant parameters. When all parameters are used (c l = 1, l = 1, . . . , 7), the WBIC value is −2324.7 and ranked at 11 (Fig. 5a-d).
WBIC is lowest (WBIC = −2327.4) when c 1 , c 2 , and c 5 = 1, which is the same combination used in the input observation series. This demonstrates that the proposed method can identify the active reaction pathways from redundant candidates of reaction pathways. For the combinations of used model parameters that minimize WBIC values (i.e., k 1 , k 2 , and k 5 ), another MCMC run was conducted to estimate the model parameters. The first 10 3 trials of MCMC sampling were not used. The auto correlation (AC) between MCMC sampling and lag time showed steep decay with increasing lag length, suggesting that MCMC sampling is less correlated and independent (Fig. 6a-c). The posterior distribution of these parameters after 10 4 Monte Carlo steps is similar to the shape of a normal distribution with a single peak (Fig. 6). The mean values of the posterior distribution for parameters log 10 (k 1 ), log 10 (k 2 ), and log 10 (k 5 ) are −3.299, −3.301, and −3.291, respectively (Fig. 6a-c), which are close to their true values.

Dependence of model selection accuracy on spatial resolution and observation noise
Here we investigated how robustly the true model can be selected for different numbers of total observations and different levels of observation noise. With a set of n and σ , the validations are repeatedly conducted for 10 trials. Different datasets are used for each trial. Figure 7a-i show the result of model selection for all ten trials at each set of n and σ . Combination of indicators ranked within 5th are shown. At n = 300, true models are definitely selected as the best (rank 1) models four times when noise levels are low (σ = 10 −2.0 ; Fig. 7a). The number of times the true models were selected are decreased at high noise levels, and are five and three times among ten trials at σ = 10 −1.5 (Fig. 7b) and 10 −1.0 (Fig. 7c), respectively. Similarly, at n = 100, the number of times the true models were selected as the best model are 2, 3, and 2 times among 10 trials for σ = 10 −2.0 (Fig. 7d), 10 −1.5 (Fig. 7e), and 10 −1.0 (Fig. 7f), respectively; At n = 30, the number of times the true models were selected as the best model are 2, 1, and 0 times among 10 trials for σ = 10 −2.0 (Fig. 7g), 10 −1.5 (Fig. 7h), and 10 −1.0 (Fig. 7i), respectively.

Dependence of parameter estimation accuracy on spatial resolution and observation noise
Here we investigated how robustly the non-zero variables (k 1 , k 2 , and k 5 ) can be estimated for different numbers of total observations and different levels of observation noise. We evaluated the discrepancy between estimated value and true value as the normalized root mean square error (NRMSE) by where k (est) l and k (true) l are the posterior mean and true value of the rate constant for reaction pathway l, respectively. For each setting of the number of observations (n) and observation noise (σ ), we evaluated the NRMSE values for several trials because the NRMSE values can be affected by a random number. Then, we used one of the trials as a typical result because the calculated NRMSE results for the several trials were consistent. Figure 8 shows that the calculated NRMSE on n and σ . We obtained NRMSE in the range of n = 15 − 300 (log 10 (n) = 1.17 − 2.48) and σ = 10 −2.0 − 10 −1.0 (corresponding to 2 − 20% deviation relative to the maximum value of each observation series). Results presented in Fig. 8 suggest that the NRMSE increases with (1) decreasing total number of observation and (2) increasing noise intensity. The NRMSE shown in Fig. 6 (log 10 (n) = 2.48 (n = 300) and σ = 10 −1.5 ) was 10 −1.7 , which corresponds to 2% error on average. Even when the number of observation was limited (log 10 (n) = 1.30 (n = 20)), rate constants were estimated with similar NRMSE (10 −1.7 ) for moderate noise levels (σ = 10 −1.75 ), whereas it was estimated with slightly higher NRMSE (10 −1.2 ) for intense noise levels (σ = 10 −1.25 ). For most regions except for log 10 (n) < 1.30 (n < 20) and σ > 10 −1.50 , the NRMSE was lower than 10 −1.0 , which corresponds to 10% estimation error on average (Fig. 8). The estimates were very accurate from the viewpoint of heterogeneous reaction in the Earth sciences, suggesting that our proposed method is effective against both observational sparseness and noise intensities and can be used for estimating kinetic parameters.

Discussion and summary
In this study, we developed a method for identifying active reaction pathways from redundant reaction pathway candidates using observed spatial data of the solid phase and Bayesian modeling. Verification using simulated data showed that the proposed method provides effective estimation for a given noisy dataset.
In general, it is difficult to determine the reaction pathways in nonlinear parallel heterogeneous reactions without measuring the solution chemistry. Even when the solution chemistry dataset is provided, the exploration and identification of reaction pathways involves multiple considerations due to the complexity and nonlinearity of the reactions. As demonstrated during the estimation of known kinetic parameters from the noisy spatial distribution of minerals, our method objectively identifies both the active reaction pathways and unimportant reaction pathways, instead of arbitrarily ignoring them. The proposed method can also be employed to evaluate the estimation accuracy for each unknown variable, which is important for scientific research.
In this study, the lowest WBIC value was not attained using all seven variables (-2324.7, corresponding to a ranking of 11; Fig. 5a). This suggests that overfitting can be avoided thorough Bayesian inversion, and the estimated model parameters are not affected by noise.
For variable combinations ranked at < 15, the indicator suggests that non-zero parameters k 1 , k 2 , and k 5 always appeared in the combinations (Fig. 5a, b); however, the number of unimportant parameters (i.e., with a zero value) included in each combinations differed. For the model ranked first, no zero-value parameters were included. For the model ranked second, the indicator suggests that k 1 , k 2 , k 5 , and k 6 were included, whereas k 1 , k 2 , k 3 , k 5 , and k 6 were included in the model ranked third (Fig. 5a, b). Because unimportant parameters were included in the models ranked second and third, the WBIC values were larger than those for the parameter combination ranked first.
Large variability was observed in WBIC values, especially between rank 48 (WBIC = 138.2) and 49 (WBIC = 5665.7; Fig. 5e). The indicators suggest that the former used three variables: k 2 , k 6 , and k 7 (Fig. 5a), whereas the latter used two variables: k 1 and k 2 (Fig.  5). When the variable combinations ranked 48 and 49 were compared, the number of variables was larger for the lower-ranked model. This suggests that the extent of fitting led to large differences in the WBIC values between the models ranked at 48 and 49.
It is important to note the characteristics of the 32 worst models identified by highest WBIC value (Fig.  5a-e). Among these worst models, all of the models without any parameters (i.e., all rate parameters set to zero) and with less than five parameter are similarly selected as worst models. This is because reactions that are necessary for another reaction to take place are restricted by these models. Common features in these worst models is that k 1 and k 2 are not used (Fig. 5a). Without these reaction pathways, the combination of variables become worst case, even when the important reaction k 5 is used. These results suggests that reac-tion pathway k 1 and k 2 are necessary to drive another reaction, such as k 5 . We suggest that identifying the common features of variables in the worst models may also be helpful to understand the sequences of reactions.
Although it is difficult to analytically derive Bayesian free energy F in general, there are several methods to calculate F [16], such as thermodynamic integration [43], nested sampling [44], and the non-equilibrium Monte Carlo method [45]. In this study, we used an information criterion, WBIC, that may provide a computationally efficient approximation of the F. The F can be approximated by a Bayesian information criteria (BIC) [31] only when the statistical model is regular, i.e., the posterior distribution can be approximated as the normal distribution [31,32]. In the example described in this study, the shape of the posterior distributions obtained by MCMC method are similar to those of the normal distribution (Fig. 6). To investigate the effectiveness of BIC, the same validation numerical experiments were conducted. As a result, we found that the BIC also identified active reaction pathways in the case of the validation test shown in the present study. However, the usage of BIC cannot be validated without sampling because the shape of the posterior distribution cannot be constrained prior to sampling. In contrast, the WBIC, which is a generalized version of the BIC, can be used in both regular and singular statistical models and regardless of the shape of the posterior distribution [32]. Analyzing real heterogeneous reactions could require a solution for a singular model because heterogeneous reactions between minerals and fluids can be complex. Thus, WBIC as a generalized version of BIC would be the preferred information criterion for selecting models via Bayesian inference.
In this study, the noise variances σ j for each observation series are treated as known constants. In the realistic situation, noise variance can be estimated as a hyperparameter by Bayesian estimation. However, such hyperparameter estimation is not necessary because approximate values for noise variance are known for the heterogeneous reaction between minerals and fluid.
In this study, spatial data points at a single time are used as observable data (n). As we showed in Fig. 7, the true model is likely difficult to be selected when n = 30 with a corresponding space interval of 0.3 cm (Fig. 7g-i). Because the spatial observation step can be 0.01 cm at minimum, n may not be less than 100, and an extreme situation in which only 30 observation data points can be used is unrealistic. However, even in the extreme situation, the parameters used (c 1 , c 2 , and c 5 = 1) in the true model frequently appeared in the top five models, suggesting that the true models can be interpreted by checking common indicators in the high ranked models.
The Bayesian approach has been applied to infer kinetic reaction pathways in the fields of chemical engineering and biochemistry [26,29]. In these works, the Bayes factor, which quantifies the support for one model over another, was used to select models. The Bayes factor between two hypothesized models is obtained by calculating the difference of F between the two models.
Such a calculation of F can be substituted by WBIC to reduce computational cost.
Although the effectiveness of our approach was demonstrated by analyzing a synthesized dataset, the limits of its applicability must be noted. It is easy to imagine that the exhaustive search method used in this study will become intractable for a large number of variables (L). Therefore, to reduce the computational load, a relaxation approach such as the least absolute shrinkage and selection operator (LASSO) method [21,46] using an l 1 -norm regularization term could be used when the ES method is computationally infeasible. However, exhaustive search (ES) is necessary for the strict selection of efficient variables [33,34,47]. This computational explosion problem associated with the ES method can be improved through the development of another relaxation method.
The proposed approach serves as a basic inversion framework for extracting active reaction pathways during parallel heterogeneous reactions. The applications of our methods to real data will be considered separately. The development of analytical methods for heterogeneous reactions and applications to a real dataset will provide a better understanding of important and dynamic Earth processes.

Author contributions
RXO conducted design of study, experiments, and drafting the manuscript. TK and TO contributed to the research idea development and revised the manuscript critically for important intellectual content. All authors gave their final approval of the manuscript version to be submitted the authors' original work, has not received prior publication and is not under consideration for publication elsewhere.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data are available from the corresponding authors upon request.] 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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Discretization of partial differential equation
The partial differential equation (Eq. 1) must be discretized to solve. The left side of Eq. (1) can be discretized with respect to time to obtain the difference equation where t denotes a time step and Δt is a time interval used for time discretization. The first term on right side of Eq.