DRGEP-based mechanism reduction strategies: graph search algorithms and skeletal primary reference fuel mechanisms

Skeletal mechanism reduction based on the directed relation graph with error propagation (DRGEP) is explored with focus on the selection of graph search algorithm and application on primary reference fuels (PRF). It is demonstrated that R-value-based breadth-ﬁrst search (BFS) and Dijkstra’s algorithm generate results independent of species order, while basic and modiﬁed depth-ﬁrst search and basic BFS depend on the order of species in the reaction mechanism. Each of the search algorithms is also used with DRGEP to generate skeletal mechanisms for a detailed mechanism for n -heptane, consisting of 561 species and 2539 reactions, covering a comprehensive range of conditions for temperature, pressure, and equivalence ratio. Dijkstra’s algorithm combined with a coeﬃcient scaling approach generates the most compact skeletal mechanism compared to the other search algorithms for a 30% error limit, with 131 species and 651 reactions. A comprehensive reduction scheme including DRGEP with Dijkstra’s algorithm followed by sensitivity analysis and unimportant reaction elimination is then applied on a large binary ( n -heptane and iso -octane) PRF mechanism, containing 1034 species and 4236 reactions, to generate skeletal mechanisms at varying levels of detail. The skeletal mechanisms with ∼ 100 species adequately predict the ignition delay and laminar ﬂame speed calculations of the detailed mechanism.


I. Introduction
T he primary obstacle to inclusion of detailed chemistry into full-scale simulations is the large and everincreasing size of kinetic mechanisms. A typical 3D internal combustion engine simulation including models for spray, combustion, and pollution is limited to about 10 radical and intermediate species (in addition to reactants and major products), 1 though this limitation is relaxed as computational hardware becomes more powerful. Surrogate models of liquid fuels of interest to transportation, aerospace, and energy applications contain blends of large hydrocarbons. However, detailed reaction mechanisms for surrogates of gasoline, 2, 3 diesel, 3,4 and jet fuels 5,6 typically include hundreds to thousands of species and thousands of reactions. For instance, a detailed mechanism for C 8 -C 16 n-alkane hydrocarbons consists of 2115 species and 8157 reactions, 7 while a recent biodiesel surrogate mechanism based on a blend of methyl decanoate, n-heptane, and methyl-9-decenoate that contains 3299 species and 10806 reactions. 8 The computational cost of chemistry scales by the third power of the number of species in the worst case when factorizing the Jacobian, 9 so models of this size are too large to include in full-scale CFD models even without considering other complex models such as spray and turbulence. In addition, the wide range of time scales associated with the large number of species and reactions induces stiffness in the governing differential equations. Reduction of the reaction mechanisms is therefore necessary to incorporate detailed chemistry into full-scale simulations.
More recently, mechanism reduction research has focused on approaches that simplify larger, detailed reaction mechanisms. Skeletal reduction is the typical first step in a comprehensive reduction scheme, where unimportant species and reactions are removed from the detailed mechanism leaving a standard kinetic mechanism. This is contrasted with time-scale reduction, which uses approximations such as quasi-steady state species and partial equilibrium reactions to reduce stiffness and requires modification to the governing differential equations. The implementation of skeletal mechanisms into CFD simulations is simpler and as such skeletal reduction is the basis of the current effort. Computational cost minimization (CCM) strategies seek to optimize operations in such a manner that introduces little or no error 9 and can be integrated with most reduction schemes. A prominent example is the in situ adaptive tabulation (ISAT) approach, 10 which stores chemical kinetics information at early stages of the simulation then later retrieves the data when needed again in order to avoid recalculation.
Though many skeletal reduction methods have been demonstrated, 9,[11][12][13] much attention recently focused on the directed relation graph (DRG) method, developed by Lu and Law. 14,15 This approach maps the coupling of species onto a weighted directed graph and performs a graph search initiating at selected target species combined with an error threshold to find unimportant species for removal. Research on DRG-based methods recently split into two directions: (1) DRG-aided sensitivity analysis (DRGASA) by Lu and Law 16,17 that follows application of DRG with sensitivity analysis to find further unimportant species and (2) DRG with error propagation (DRGEP) by Pepiot-Desjardins and Pitsch 18 that considers the propagation of species dependence down graph pathways. Both modifications showed noticeable improvement in the magnitude of reduction over the original DRG method, though DRGEP has seen wider usage. [18][19][20][21][22][23][24] Most recently, Niemeyer et al. 22 presented a combination of DRGEP and DRGASA called DRG with error propagation and sensitivity analysis (DRGEPSA) that showed greater improvement in reduction over all three preceding methods. Sun et al. developed another DRG-based method, path flux analysis, 25 that uses both production and consumption fluxes to define the directed graph and identify important species but is unable to achieve the level of reduction of methods that include sensitivity analysis. All the above-mentioned DRG-based approaches demonstrated the potential of reducing large detailed mechanisms for higher-order hydrocarbon fuels into more manageable skeletal mechanisms.
The DRG method maps the coupling of species onto a directed graph and finds unimportant species for removal by eliminating weak graph edges using an error threshold, where the graph nodes represent species and graph edge weights indicate the dependence of one species on another. Following a simple graph search initiated at certain preselected target species (e.g., fuel, oxidizer, important pollutants), species unreached by the search are considered unimportant and removed from the skeletal mechanism. DRGEP modifies this approach slightly by considering the dependence of a target on other species due to a certain path as the product of intermediate edge weights and the overall dependence is the maximum of all path dependencies.
One point that has not received much attention is the choice of graph search algorithm used to calculate this overall dependence. While the DRG method only needs to find connected graph nodes and can therefore use any search algorithm, caution must be taken when selecting the method used with DRGEP. The issue of calculating the overall dependence is actually a single-source shortest path problem 26 where the "distance" of the path is the product of intermediate edge weights rather than the sum, and the "shortest" path is that which has the maximum overall dependence rather than the minimum. Search algorithms that in general do not correctly calculate the overall dependence will underestimate the importance of species and cause premature removal from the skeletal mechanism. Further, the resulting skeletal mechanism is independent of the order of the species in the mechanism only if a specific algorithm can correctly capture and calculate the overall dependence of species. Therefore, the reliability of various algorithms needs to be studied by determining whether each is dependent on the order of species in a detailed mechanism.
Most DRGEP studies reported using algorithms based on either depth-first search (DFS) 22 or breadthfirst search (BFS), 20,21,23,24 but no comparison has been performed. Here we compare such methods with Dijkstra's algorithm, the classical solution to the single-source shortest path problem, 26,27 in order to demonstrate the weaknesses of DFS-and BFS-based approaches and the subsequent effectiveness and reliability of Dijkstra's algorithm. First, by randomly shuffling the order of species in the detailed mechanism for n-heptane of Curran et al., 28 we show that DFS-and BFS-based algorithms can generate results dependent on the order of species in the reaction mechanism while Dijkstra's algorithm generates consistent results regardless of species order. Second, we demonstrate that, for a given error limit, more compact skeletal mechanisms are possible with DRGEP by using Dijkstra's algorithm. This is done by comparing skeletal mechanisms generated with different search algorithms for n-heptane covering a comprehensive range of au-toignition conditions for pressure, temperature, and equivalence ratio. Third, we demonstrate the reduction capability of DRGEPSA using Dijkstra's algorithm on multicomponent fuel blends by generating compact skeletal mechanisms from the detailed n-heptane/iso-octane Primary Reference Fuel (PRF) mechanism of Curran et al. 29,30 containing 1034 species and 4236 reactions.

II.A. DRGEP Method
The skeletal reduction approaches used here are based on DRGEP and DRGEPSA, described in detail by Pepiot-Desjardins and Pitsch 18 and Niemeyer et al., 22 respectively. Accurate calculation of the production of a species A that is strongly dependent on another species B requires the presence of species B in the reaction mechanism. This dependence is expressed with the direct interaction coefficient (DIC): where A and B represent the species of interest (with dependency in the A → B direction meaning that A depends on B ), i the i th reaction, ν A,i the stoichiometric coefficient of species A in the i th reaction, ω i the overall reaction rate of the i th reaction, and n R the total number of reactions. After calculating the DIC for all species pairs, a graph search is performed starting at user-selected target species to find the dependency paths for all species from the targets. A path-dependent interaction coefficient (PIC) represents the error propagation through a certain path and is defined as the product of intermediate DICs between the target T and species B through pathway p: where n is the number of species between T and B in pathway p and S j is a placeholder for the intermediate species j starting at T and ending at B. An overall interaction coefficient (OIC) is then defined as the maximum of all PICs between the target and each species of interest: Pepiot-Desjardins and Pitsch 18 also proposed a coefficient scaling procedure to better relate the OICs from different points in the reaction system evolution that we adopt here. A pseudo-production rate of a chemical element a based on the production rates of species containing a is defined as where N a,S is the number of atoms a in species S and P S and C S are the production and consumption rates of species S as given by Eqs. 2 and 3, respectively. The scaling coefficient for element a and target species T at time t is defined as For the set of elements {E}, the global normalized scaling coefficient for target T at time t is Given a set of kinetics data {D} and target species {T }, the overall importance of species S to the target species set is We employ the same sampling method as given by Niemeyer et al., 22 where constant volume autoignition simulations are performed using SENKIN 31 with CHEMKIN-III. 32 Species are considered unimportant and removed if their R S value falls below a cutoff threshold ε EP , which is selected using an iterative procedure based on a user-defined error limit for ignition delay prediction. 22 Reactions are eliminated if any participating species are removed.

II.B. Graph search algorithms
In this study we compare the results of DRGEP using basic DFS, modified DFS, basic BFS, R-value-based BFS (RBFS), and Dijkstra's algorithm. Cormen et al. 26 presented detailed discussion of and pseudocode for DFS, BFS, and Dijkstra's algorithm, while the modified DFS used by Niemeyer et al. 22 and RBFS of Liang et al. 20 differ only slightly from the basic DFS and BFS.

II.B.1. DFS-based algorithms
The DFS initiates at a root node, in this case a target species, and explores the graph edges connecting to other nodes. The first node found is added to a last in, first out stack, then the search moves to this node and repeats. In this manner, the search continues deeper down the graph pathway until it either reaches a node with no connections or all the connecting nodes have been explored, then backtracks one position up the stack. The search is performed separately using each target species as the root node and the maximum OIC is stored for each species. Lu and Law first introduced the DFS in this context for the DRG method. 14 The modified DFS used by Niemeyer et al. 22 for the DRGEP method (but not described in detail) differs from the basic DFS in that the OIC values for all target species are set to unity before starting the search at the first target. The search then only repeats starting at other target species not discovered in the initial search. The resulting OIC values combine the dependence of all targets on the species and this prevents the use of coefficient scaling based on individual target species activity.

II.B.2. BFS-based algorithms
The BFS initiates at the root node and explores all adjacent graph edges, adding connected nodes to a first in, first out queue. After discovering all connected nodes, the search moves to the first node in the queue and restarts the procedure, moving to the next node in the queue after discovering all connected nodes. Previously discovered nodes are not repeated. As with the DFS, the search initializes at each target separately and the maximum OIC for each species is stored.
The R-value-based BFS (RBFS), first introduced by Liang et al. 20,21 and also used by Shi et al., 23,24 differs from the standard BFS in that PICs smaller than the error threshold ε EP are not explored. In other words, the search ignores graph pathways that would lead to an OIC below the cutoff value and only allows discovery of important pathways. This increases efficiency by avoiding exploration of unimportant pathways but causes the RBFS to depend on the value of ε EP , unlike the other search methods. As such, the primary application of this algorithm lies in dynamic reduction where the threshold value is known a priori, rather than general comprehensive skeletal reduction where the threshold is to be determined iteratively based on a user-defined error limit. 22

II.B.3. Dijkstra's algorithm
Dijkstra's algorithm, originally introduced by Dijkstra 27 and discussed in detail by Cormen et al., 26 differs from DFS, BFS, and their variants because it was specifically designed to calculate the shortest paths from a single source node to all other nodes rather than simply search the graph to find connected nodes. As mentioned previously, calculating OIC values is a shortest path problem where the "shortest" path is that with the maximum product of intermediate edge weights. This is the only modification needed to apply Dijkstra's algorithm to the current problem. In brief, Dijkstra's algorithm functions similarly to BFS except that it starts with the set of graph nodes stored in a max-priority queue. The algorithm finds and removes the node with the maximum OIC value (initially, the root node) from the queue and explores adjacent nodes to calculate or update OIC values. When all neighboring nodes are explored, the procedure restarts at the node with the next-highest OIC in the queue. The search completes when the queue is empty. As with previously-discussed search algorithms, this is performed separately for each target species as root node.

II.C. Sensitivity analysis
In DRGEPSA, sensitivity analysis (SA) is performed by first comparing the OIC values of remaining species with an upper threshold ε * (e.g. 0.2-0.4). Species where ε EP ≤ R S < ε * are labelled as "limbo" to be further analyzed for removal, while species where R S ≥ ε * are automatically retained in the final skeletal mechanism. Limbo species are removed from the mechanism one-by-one then sorted using: where δ B,ind is the induced error due to the removal of species B with respect to the detailed mechanism and δ DRGEP is the error of the DRGEP mechanism. The limbo species are then removed in order until the error reaches the user-specified error limit.

II.D. Unimportant reaction elimination
Following the DRGEP and SA phases, an additional step of further unimportant reaction elimination is performed based on the methodology of Lu and Law. 17 Using an approach based on the CSP importance index, 33 the normalized contribution of a reaction i to the production rate of a species A is where a reversible reaction must be treated as a single reaction. 17 reaction i is retained in the mechanism. The threshold ε reac is determined iteratively based on a user-defined error limit in a similar manner to ε EP .

III.A. Reliability of graph search algorithms
In order to demonstrate the dependence of DFS-and BFS-based methods on the order of species in a reaction mechanism, all five graph search algorithms are used with DRGEP to generate skeletal mechanisms from the detailed mechanism for n-heptane of Curran et al., 28 containing 561 species and 2539 reactions, where the order of species is randomly shuffled. Basic DFS, basic BFS, RBFS, and Dijkstra's algorithm are used both with and without the target-based coefficient scaling while the modified DFS is not compatible with this and so is used without it. Autoignition chemical kinetics data are sampled from initial conditions at 1000 K, 1 atm, and equivalence ratios of 0.5-1.5. Oxygen, nitrogen, n-heptane, and the hydrogen radical are selected as target species. Figure 1 shows the ignition delay error of skeletal mechanisms at varying levels of detail generated by DRGEP using RBFS and Dijkstra's algorithm with coefficient scaling and the modified DFS. Basic DFS and BFS demonstrate similar dependence on species order to the modified DFS and thus are omitted from Fig.  1 for clarity, but the skeletal mechanism sizes from each may be different as will be shown. Similarly, the coefficient scaling comparison is also omitted because its application does not affect dependence on species order, although the resulting skeletal mechanism sizes are different. Comparison of results from the original mechanism and five mechanisms with randomly shuffled species illustrates the dependence of the modified DFS on the order of species in the mechanism while RBFS and Dijkstra's algorithm produce consistent results regardless of species order. This is because DFS and BFS explore graphs based on the order of nodes while Dijkstra's algorithm follows the strongest path independent of order. 26 In addition, Fig. 1 shows that both Dijkstra's algorithm and RBFS produce smaller skeletal mechanisms before the error becomes unacceptably large. The dependence of DFS, modified DFS, and BFS on species order stems from the fact that these algorithms do not in general calculate the correct OIC for every species and as such can underestimate the importance of species, causing unwarranted removal. RBFS produces the same results as Dijkstra's algorithm because it prevents exploration of unimportant paths, i.e. paths with PICs (see Eq. 5) less than the threshold, and therefore only discovers paths that lead to an OIC greater than the threshold value. While this value may be smaller than the correct OIC as calculated by Dijkstra's algorithm, the species is nonetheless considered important and retained in the skeletal mechanism. This subtle error could be significant, though, if the DRGEP method is followed by a further reduction stage such as sensitivity analysis that depends on the OIC values. 22

III.B. Effectiveness of graph search algorithms
Skeletal mechanisms for n-heptane are generated covering a comprehensive range of conditions using basic DFS, modified DFS, basic BFS, and Dijkstra's algorithm. Basic DFS, BFS, and Dijkstra's algorithm are used with and without the coefficient scaling in order to determine its effect as well. Autoignition chemical kinetics data are sampled from initial conditions at 600-1600 K, 1-20 atm, and equivalence ratios of 0.5-1.5. Oxygen, nitrogen, n-heptane, and the hydrogen radical are selected as target species, and the error limit in ignition delay prediction is 30%. RBFS is not compared here because it is not suited for a comprehensive skeletal reduction due to its dependence on threshold value; additionally, based on the results in the previous section we assume the resulting skeletal mechanism would match that from Dijkstra's algorithm.
The results using the original mechanism of Curran et al. 28 are summarized in Table 1. All search methods lead to skeletal mechanisms with similar performance, but Dijkstra's algorithm with coefficient scaling produces the most compact skeletal mechanism. BFS and modified DFS produce similar results while basic DFS is unable to generate a comparable skeletal mechanism for the given error limit. This weakness is most likely due to the fact that basic DFS finds only the first path from target to species, which typically contains many intermediate species and severely underestimates the OICs for some important species. The modified DFS performs better by inserting the other targets into the pathways and therefore increasing the OIC values for important species, while BFS finds the path with the shortest number of intermediate nodes. 26 Table 1 also demonstrates the need for coefficient scaling. Both basic DFS and Dijkstra's algorithm generate more compact skeletal mechanisms with the scaling, while BFS actually produces a slightly larger skeletal mechanism. Without scaling, Dijkstra's algorithm is unable to produce a more compact skeletal mechanism than the other methods, despite its ability to generate results independent of species order in the parent mechanism.

III.C. PRF skeletal mechanisms
Skeletal mechanisms at varying levels of detail are generated from the detailed mechanism for binary PRF mixtures (n-heptane and iso-octane) of Curran et al., 29,30 containing 1034 species and 4236 reactions, using DRGEPSA (with Dijkstra's algorithm) followed by unimportant reaction elimination. Two sets of constant volume autoignition initial conditions are used to generate chemical kinetics data: 1) 600-1600 K, 1 atm, and equivalence ratios 0.5-1.5 and 2) 1000-1600 K, 1 atm, and equivalence ratios of 0.5-1.5; the former is referred to as the comprehensive-temperature range and the latter as the high-temperature range. The PRF mixture used has an octane number (ON) of 80, corresponding to 80% iso-octane and 20% n-heptane liquid volume at room temperature. Error limits of 10%, 20%, and 30% are used and the DRGEP target species are n-heptane, iso-octane, oxygen, and the inert nitrogen (to prevent removal). Tables 2 and 3 summarize the skeletal mechanism results for the comprehensive-and high-temperature condition ranges, respectively. As expected, in general with increasing acceptable error limit the numbers of species and reactions in the skeletal mechanisms decrease, although the 20% high-temperature skeletal  Table 2: Skeletal mechanisms at multiple levels of complexity generated from the detailed mechanism for binary PRF mixtures of Curran et al. 29,30 using DRGEPSA followed by unimportant reaction elimination. The kinetics data are generated from n-heptane/iso-octane (ON 80) autoignition at 600-1600 K, atmospheric pressure, and varying equivalence ratios.
mechanism is actually larger than the 30% skeletal mechanism. This is due to the generally nonlinear relationship between number of species and error in ignition delay prediction. The high-temperature skeletal mechanisms are more compact than the comprehensive-temperature skeletal mechanisms, consistent with previous results for n-decane. 22  Table 3: Skeletal mechanisms at multiple levels of complexity generated from the detailed mechanism for binary PRF mixtures of Curran et al. 29,30 using DRGEPSA followed by unimportant reaction elimination. The kinetics data are generated from n-heptane/iso-octane (ON 80) autoignition at 1000-1600 K, atmospheric pressure, and varying equivalence ratios.

Error limit
Validation of the skeletal mechanisms is performed with autoignition simulations over a range of initial conditions for temperature, pressure, and at varying equivalence ratios, using the PRF mixture with octane number 80 (referred to as PRF80) as well as neat n-heptane and iso-octane. The results for the comprehensive-temperature skeletal mechanisms are presented in Fig. 2 and for the high-temperature mechanisms in Fig. 3. All the skeletal mechanisms well predict ignition delay compared to the detailed mechanism, with the largest errors at high pressure (likely due to atmospheric pressure being used for the reduction kinetics data). The most noticeable error in the prediction of the comprehensive-temperature skeletal mechanisms is in the negative temperature coefficient (NTC) region, as well as at lower temperature and high pressure for n-heptane. The high-temperature skeletal mechanisms similarly display the largest discrepancy at high pressure and lower temperatures, due to the removal of low-temperature chemistry.
Further validation of the PRF skeletal mechanisms is performed with laminar flame speed calculations at   atmospheric pressure with an unburned mixture temperature of 400 K, displayed in Fig. 4. The DRGEP-stage high-temperature skeletal mechanisms are included for further analysis in Fig. 4b. All the skeletal mechanisms well reproduce the laminar flame speeds of the detailed mechanism except both the comprehensiveand high-temperature 30% error limit skeletal mechanisms, which grossly miscalculate the entire curve especially at lean equivalence ratios. The good results of the DRGEP 30% high-temperature skeletal mechanism suggest that the sensitivity analysis phase eliminated some species necessary for flame calculations but unimportant to the ignition delay. Detailed study finds that sensitivity analysis eliminated the species ic3h5co from both the comprehensive-and high-temperature skeletal mechanisms, and adding it back corrects this large error in both. Although a tighter error limit prevents the removal of this important species, this result justifies either the need for a posteriori validation or inclusion of further targets for error evaluation (i.e., laminar flame speed). Also noteworthy are the good predictions of the DRGEP high-temperature skeletal mechanisms, further validating the use of DRGEP in dynamic skeletal reduction approaches. 20,21,23 Additional laminar flame speed validation of the high-temperature skeletal mechanisms is performed by varying the ON from 0-100 at stoichiometric conditions and atmospheric pressure, with unburned mixture temperatures of 400 K and 450 K. The results for the DRGEP and final skeletal mechanisms are displayed in Fig. 5. Similar to the previous flame speed results, all the skeletal mechanisms perform well for the full ON range except the 73-species mechanism. When corrected with ic3h5co, this mechanism also demonstrates good performance similar to the others. One interesting result is that the 73-species mechanism does predict the ON 0 (n-heptane) flame speed accurately, though it fails with increasing iso-octane content.

IV. Conclusions
Graph search algorithms used in the DRGEP method for skeletal mechanism reduction are compared. Basic DFS, basic BFS, modified DFS, RBFS, and Dijkstra's algorithm are implemented in DRGEP and used to generate (1) skeletal mechanisms covering a limited range of conditions using a randomly shuffled detailed mechanism for n-heptane to determine the dependence of results on species order and (2) skeletal mechanisms covering a comprehensive range of conditions to determine the effectiveness of each algorithm. The RBFS algorithm is not used to generate a comprehensive skeletal mechanism because it is not suitable for use when the cutoff threshold is not known a priori. Both Dijkstra's algorithm and RBFS are able to generate consistent results independent of species order, and Dijkstra's algorithm used with target-based coefficient scaling generates a more compact skeletal mechanism than the other methods. Even with the improved search algorithm, however, the size of the resulting skeletal mechanism (131 species and 651 reactions) is larger than that of DRGEP with sensitivity analysis (108 species and 406 reactions), 22 suggesting that the post-DRGEP sensitivity analysis is still required.
Skeletal mechanisms for the binary PRF mechanism of Curran et al., 29,30 containing 1034 species and 4236 reactions, are also generated to demonstrate the reduction capability of a DRGEP-based comprehensive reduction scheme, where DRGEP (using Dijkstra's algorithm) is followed by sensitivity analysis and then fur- ther unimportant reaction elimination. Skeletal mechanisms covering comprehensive-and high-temperature conditions at 10%, 20%, and 30% error limits are demonstrated and validated using autoignition and laminar flame simulations. All the skeletal mechanisms demonstrate good performance in reproducing the global characteristics of the detailed mechanism, except both the 30% skeletal mechanisms which fail to reproduce the laminar flame speed especially at lean equivalence ratios. This error is attributed to the removal of the species ic3h5co by sensitivity analysis, and the addition of this species back to both mechanisms corrects it. A guided sensitivity analysis may prevent such a removal. While the comprehensive reduction scheme applied here significantly decreases the size of the large multicomponent PRF mechanism (1034 to ∼100 species), further reduction stages are necessary to integrate this detailed chemistry into full-scale CFD simulations. The problem is further exacerbated by the everincreasing size of multicomponent surrogate fuel mechanisms. Isomer lumping and diffusive species bundling are two stages currently being investigated which can offer additional computational speedup without increasing error, 17 particularly for large hydrocarbons. In addition, dynamic reduction techniques based on DRGEP similar to that of Liang et al. 20,21 and Shi et al. 23 are also being pursued, which requires further development but is potentially the most promising approach.