Node-of-Influence Network Analysis for Targeted Condition Sequencing in Plasma Chemical Reaction Networks

It was demonstrated that a plasma chemical reaction system can be represented as a directed bipartite variable-relationship (VR) graph to accurately represent node influence, and that targeted node-influence analysis of this graph can determine improved condition sequencing for a target outcome. A novel graph algorithm (OCARINA) was adapted for VR-graphs to give a measure of the net-influence of incremental increase (NIII) of one variable on another at any depth in the graph. Additionally, two conventional node-influence measures, the Eigenvector Centrality Index (ECI) and Katz Centrality Index (KCI), were also trialled on the VR-graph. The electron energy (ε) node influence was evaluated on a “baseline” continuous sequence of 10 ns 1 eV pulses in a 0D chemical-kinetic simulation using ECI, KCI and OCARINA NIII at three depths. KCI appeared to give meaningful values for ε influence in the whole graph but not on specific nodes, ECI gave no meaningful results. OCARINA (O3 targeted) suggested each successive ε pulse had diminishing influence on O3 formation, though analysis differed for each NIII depth. O3 concentrations in simulations with different numbers of ε pulses decreased with each additional pulse, correlating with the OCARINA analysis. O3, NO and O Species production in simulations of two ε pulses with one or both pulses incrementally changed by 10% from the baseline also largely agreed with the OCARINA results for each species on a baseline simulation of two consecutive 1 eV electron energy pulses. Additionally, it was found the NIII at antecedent depths corresponded to effects in subsequent phases in simulations.


Condition sequencing-The Staging of Input Conditions for a Desired Chemical Output
There is an ever-increasing variety of different plasma sources and power supply options for plasma chemistry, producing a range of different plasma discharges [1]. Though this results in a widening variety of different options, each option ultimately provides differing conditions (e.g. electric field strength, temperature, feed gas flow rate) in different sequences (e.g. 1000 °C then − 50 °C then 25 °C etc.) to an Eularian spatial element of the plasma volume.
Which sequence of conditions is most appropriate depends on the required outcome (formation of a particular chemical, for example), and if an ideal condition sequencing for a given target outcome can be identified, this can be used as a basis for selecting the most suitable options to provide them.
For electron energy, for example, if the ideal sequence of electron energy can be known, it provides a focus for the required timing of voltage amplitudes, and therefore which power supply settings to use for the plasma source, thus narrowing down the vast array of potential options to a smaller number of more suitable contenders. In addition to electron energy (which tends to be the main controllable 'input' in non-thermal atmosphericpressure plasma sources), there are also points when specific levels of chemicals, catalysts, heat, pressure, and almost any other variable in the system would be conducive to a desired outcome and could (if practical) also be sequentially altered accordingly. This is also true for almost any chemical process. For a simple process or system, finding these points is essentially trivial. In chemical processing, the staging of operations in time or space-time is easy to represent and optimise for small reaction systems, the fed batch reactor or the intermediate inlets in a tubular reactor are well-known examples of designing by this principle [2] and for reactors-in-sequence, the Levenspiel plot [3]. But as processes/systems become increasingly complex the means of analysing them to find the ideal staging/sequencing of conditions need to be increasingly sophisticated. For networks with many reaction pathways and intermediate species, it is an open question. Presented here is an attempt to address it for plasma chemical reaction networks.
Increasingly complex systems mean increasingly greater numbers of entities which can influence each other in an increasing number of ways. Fortunately, it is possible for these influences between entities to effectively be mapped as a mathematical graph/network (the terms graph and network are used interchangeably in this work), and the resulting network topology to be analysed.
Mathematical graphs are composed of nodes and edges, with nodes being able to represent a variety of objects and edges representing the relationship between such objects. These graphs can hold a wealth of information, and additional information can be provided alongside the nodes and their relationships, to aid analysis. For example, it is possible to attribute direction [4] and weight to edges within a graph, as used in this research. Graph theory encompasses the study of these mathematical graphs and their properties, as well as algorithms that can be used to extract information from these graphs.
Graph theory has long been used in chemistry. A popular application is to model the synthesis of molecules [5,6]. A further application is to model chemical reactions using graphs (it is within this area that our own research lies), and there are many approaches to doing this, including modelling exclusively the chemical species as nodes [7][8][9][10]. Other approaches separate the nodes into two distinct categories: species and reaction [11][12][13][14] and using edges to connect reactants to a specific reaction node for any given reaction, and then connecting the product nodes to this reaction as demonstrated in Holmes et al. [15], for application to plasma chemical reaction engineering. These species-reaction graphs fall into the category of a bipartite graph since they have two distinct sets of nodes (species and reaction) such that no node in one set is adjacent (connected by an edge) to any given node in the same set, only to nodes in the other set.
Reaction networks can also be modelled as hypergraphs [16]. These are graphs such that an edge does not only have to connect two nodes at a time. An edge might exit one node and partition into n edges which then connect to n new nodes. Fagerberg et al., [6] demonstrate that the directed graphs they use to model synthesis can be represented as hypergraphs. Current software packages, however, do not adequately support the visualisation and manipulation of hypergraphs for complex graphs (such as those we demonstrate in this work), though hypergraphs may be well worth considering for solving these types of problems in pending further developments in this area.
In plasma chemistry, graph theory has been used to filter reaction sets down to reducedorder plasma chemical reaction systems that can be used in the place of larger ones to significantly reduce computation time without losing any key functionality [17]. This is extremely useful for incorporating plasma chemistry into larger models (such as 2D or 3D time-dependent plasma simulations). Whilst not explicitly employing graph theory, there are other methods that deal with chemical reaction pathways based on the methods developed by Lehman [18] and adapted specifically to plasma chemistry with PumpKin [19]. These methods have been compared to graph theoretical alternatives [20] for identifying the principal pathways relating to a given species, and eliminating negligible pathways, reactions and species for computational expediency. Whether or not to apply such reduction methods depends on the purpose of modelling the chemistry. In chemical engineering, and engineering in general, one key use of modelling a system is to find the optimal conditions for a desired outcome. There is often a trade-off between having a model simple enough to do the number of parametric sweeps necessary to find these conditions, and having a model accurate enough that the results of such parametric sweeps will apply in reality. Even models with reduced reactions sets can still be complex enough to make a complete parametric sweep computationally unfeasible, and in the process, important information about potential key reactions may be lost, especially for reduction methods which eliminate reactions based on the roles they have been previously observed to play. Two chemical species, for example, may be extremely scarce, but when in the vicinity of each other may have a high probability of interaction (large collisional cross section) and form valuable products. Since this reaction would not contribute much to one or more simulations it might therefore be eliminated from future simulations for computational expediency. But since in many chemical engineering scenarios it is possible to inject reactants at any point in a process if necessary, or to catalyse reactions, if any such reactions have the potential to significantly influence a desired outcome, the obstacle of their low concentration of reactants may seem like an obstacle worth overcoming. Similarly, if any slow reactions form products that are potential key species for a desired outcome, means of increasing their rates may be worth looking into, even if other reactions tend to predominate in conventional conditions. In graph theory, a reaction (or species) such as this, if represented as a node, is a potential "node of influence" either on another node or number of nodes, or on the graph as a whole. Reactions and species with such potential can be elucidated by making plasma species-reaction graphs with edges weighted according to the reaction rate coefficients, as suggested in [15].
A node of influence is itself a flexible term. It could mean a node that directly influences another node, a node that influences the whole graph, a node that influences a subgraph of the graph, or a node that influences another node via its influence on the whole graph. The latter is of particular interest, but there do not seem to be many algorithms available that achieve this. This may partly be due to the types of graph used. Finding nodes of influence relies on having a graph where the influence of nodes on each other is accurately represented to begin with.
This work looks at the analysis for node-of-influence on target nodes and how this influence changes over the course of a 0D chemical kinetic simulation, with the aim of assessing how to alter these nodes to find improved condition sequencing for a target outcome. First, a more appropriate graph to represent node influence will be described in the next section. Following this, some key node-influence algorithms are discussed, including a modified version of the one in [15]. The next section describes the results of an investigation where a 0D simulation was run and the results of each time step were used to update a specific graph type on which, in turn, some node influence algorithms were run at regular intervals during the course of the simulation. Following these initial results, the results of the modified version of the algorithm in [15] were compared to the results of 0D simulations of different condition sequences (specifically of electron energy peaks) to assess the effect of altering nodes of influence according to the results of this node-influence algorithm.

Variable-Relationship graph-An alternative graph format to more accurately represent influence
Previously Species-Reaction graphs were used in Holmes et al. [15]. However, a speciesreaction graph falls short of an accurate representation of all the nodes of influence for plasma chemistry, as arguably one of the greatest sources of influence in a plasma chemical system is the electron energy/electron temperature, which is neither a species or a reaction, and lies outside a species-reaction graph. The method described in the following section is an advance on this. In a simple mathematical description of chemical reactions (often used in reaction engineering), the species concentrations are variables, and the relationships between these variables are the reaction rate equations. Therefore, bipartite graphs of species and reactions can be generalised to graphs of variables and relationships. In the following example, for a reaction R which is the rate coefficient (k) is often a function of the electron temperature (T e ) in a plasma chemical reaction set, and the relationship of the change in concentration for C and D would be described by the following equation: And the relationship of the change in concentration for A and B would therefore be: A graph can be made from this relationship, with electron temperature also added as a variable node in the graph. This is shown in Fig. 1.
The node R in the centre of Fig. 1 is a relationship node, all other nodes in the figure are variable nodes and hence the graph is a variable-relationship graph (VR graph). Whether or not a variable is a dependent variable or an independent variable with respect to the relationship is represented by the direction of the edge (VR is an independent variable, and RV is a dependent variable). Due to the nature of equations expressing the concentration changes of reactants, it is necessary for the graph to be allowed negative edge weights as well as positive ones (e.g. edges RA and RB in Fig. 1).
The edge weights should represent the influence of the start node on the end node. For edge (A, R) or AR in Fig. 1, for example, the edge should constitute some form of description of the influence variable A has on relationship R. Clearly, the only effect of any influence on a scalar quantity is to increase or decrease its value, and the value of an equation will change depending on the value of its independent variables, and this value in-turn determines the value of the dependent variables, which are also independent variables of successive relationships, and so on.
A variable-relationship graph can be evaluated from the equations of its relationships to represent the system in one or more states where all the variables have specific values, or alternatively it could represent some measure of the range of values (e.g. minimum to maximum) which each variable and relationship can take. It can also be used without being evaluated simply to investigate the connections between the variables and relationships irrespective of their values (for simple measures of node connectivity, for example).
For graphs where the variables have specific values, edge weights can be set according to the partial differential of the start node with respect to the end node (Fig. 2). The result is a dynamic graph with the edge weights representing the instantaneous influence each node has on its successor. In the case of relationship-to-variable edges, the infinitesimal change in the dependent variable is identical to the infinitesimal change in the relationship, R, or identical to -R. Thus the ratio of the change in the value of the relationship to the change in the value of the dependent variable, and thus the edge weight, Another key piece of information is whether the influence of a variable on a relationship is always positive or always negative over a given value range, or always high or always low in relation to a defined value over a given range. Any of these four cases effectively means that although the value of the variable and relationship may change, this specific property (e.g. positive or negative, relatively high value or relatively low value, etc.) of the V → R edge in the system does not change for any possible scenario. Each of the plots in Fig. 4 gives some examples of range values. Such range graphs, can be used for algorithms where the system can be accurately represented over the given ranges. It may be necessary to use more than one number to accurately summarize the relationships in Fig. 4, (especially as the relationships are frequently functions of more than one variable), but various graph libraries (e.g. networkx) allow multiple edge attributes to be specified, so this should not necessarily be prohibitive.
It may also sometimes be useful to visualise these networks. In a network/graph visualisation, there are a number of scales that each of these values could be assigned to, e.g. Edge thickness, Edge label, Edge label size, Colour Hue, Colour saturation, Colour brightness. This can be done using software such as Gephi [21] as demonstrated in [15].  In addition to the tools available for visualisation, there is a wealth of algorithms available that could be used to analyse directed graphs and most programming languages have dedicated libraries for this purpose. Such algorithms are used to analyse mathematical graphs when the solution cannot be directly calculated, and as such, require heuristic methods or approximation algorithms. If such problems are hard to solve, but easy to verify a potential solution once it has been found, the problem is said to be NP-complete. The problem of maximising output in a chemical reaction network has been shown to be an NPcomplete problem [22]. Regarding such problems, Fagerberg et al., [6] make the point that several good plans to choose from are better than one. Further to this, if all the options for improving a chemical synthesis can be presented (e.g. an algorithm could present them), then improvements to a current or baseline process can be made according to one's judgement and the resources one has available, to a much greater extent than simply selecting a plan. In the context of this work, a plan can be thought of as an assembly of elements made of expected causes and effects (X causes Y, Y causes Z, etc.). If it is possible to consider all the elements of cause and effect then all the options of what could be conducive or destructive to a single objective could be laid out, allowing options to be selected that are best suited to available resources/capabilities. A graph can serve this purpose. The challenge is in identifying which elements are important to include in this "palette" of options, and which are not, and how to include them. Visualisation can be an important component in the ability to do this, but cannot be solely relied upon, especially when graphs become increasingly complex.

Graph analysis methods for nodes-of-influence
There are many existing graph algorithms which are all variations of centrality measures which rank nodes according to a measure of their importance within the graph. There are a variety of such measures, and they work to varying degrees of success depending on the graph they are being used on. Two of the most widely used of these methods are the Eigenvalue centrality index, and the Katz centrality index, which are both measures of node centrality that consider the importance of neighbour nodes as well as the node itself. The importance of a node increases if it is itself connected to an important node [23]. This could be particularly convenient given the sequential nature of chemical reactions. For any given graph, a matrix A can be formed such that entry a i,j in the matrix is equal to one if there is an edge connecting node i and j . This matrix is called an adjacency matrix and can be further extended to weighted graphs by replacing the value of one with the weight of the edge connecting i and j . The adjacency matrix can be used in measures of centrality. In graph theory, measures of centrality allow us to determine how influential any given node is, relative to other nodes in the graph. Measures of centrality have been used for highlighting a node's importance in chemical reactions, such as Silva et al. [14], who use the outdegree centrality to highlight influential nodes. However, this measure does not take into account if a node is connected to other 'influential' nodes and so we argue that eigenvector centrality measures make for a more natural metric to use.
By taking the eigenvector that corresponds to the highest eigenvalue of the adjacency matrix, you obtain a weighted sum of not only the direct connections but the indirect connections of any given node [23].
For a graph G ∶= (V, E) with adjacency matrix A = a i,j the relative centrality score x i of node i is given as: such that M(i) is the set of neighbours for node i and is constant. There are multiple variations of this measure. One such variation is the Katz Centrality index. In this case, the Katz centrality measure [24] is given by: such that and are positive constants and contributes to all vertices and removes the possibility of a node being attributed null centrality if it is not strongly connected to other components of the graph.
Another widely used centrality method is the page rank algorithm, which is a subtle adjustment to the Katz centrality. This, along with the eigenvector and Katz centralities are all described in [25].
These centrality measures, however, are concerned with node influences on the whole network. If we are seeking to target a particular chemical species represented as a node, we are interested in node influences on a specified node via the whole network. Figure 5 outlines a variable-relationship graph of a simple hypothetical reaction network composed of 3 chemical reactions: R 1 = A B + C, R 2 = B + C E and A + B A + D, and how this graph would be operated on using this process. In the example in Fig. 5, D is the target node. How important each node is to the target node depends on how it positively or negatively influences the target node, either directly (if it is a predecessor) or indirectly (if it is a predecessor of a predecessor… etc.). Using Variable-Relationship graphs of the type specified in Fig. 2, this importance can be determined numerically.
This has been investigated to some degree in [15]. Ranged species-reaction graphs were also used in a graph algorithm aiming to make the first few steps towards a graph theory based waveform targeting algorithm based on the ranges of the rate coefficients over a given range of electron energies, which the authors named Optimal Condition Approaching via Reaction-In-Network Analysis (OCARINA). This was based on the following logic: "If a reaction creates the target it's desirable, if a reaction destroys a target it's undesirable." And by extension: "If a reaction creates something that creates the target it's desirable, and if a reaction destroys something that creates the target, it's undesirable, etc.". The inclusion of negative edges (see Figs. 1, 2 and 3) in the VR graph, however, means this distinction of "negative and positive" is already made by the graph edges. Thus the influence of any node (the predecessor) on any other node (the target) can be measured by the sum of all the edges connecting the predecessor to the target, which can be performed on the adjacency matrix. The result will be a vector of the node influences (node ranking vector). The adjacency matrix can be multiplied by this vector, and summed over again for successive "depths" of predecessor any number of times ( Fig. 5), which can be summarized as where x i is the ith column of the node ranking vector, a i,j is the jth row of the ith column, and n is the depth of predecessor from the target node. The initial node ranking vector (n = 0) with the index of each node set to 0, except the target node which is set to 1. "Depth" refers to the number of "steps back" (equal to the number of edges) between a predecessor node and the target node [26]. This method is explained in much more detail in the Appendix, including a validation of how the if-then statements in [15] determining the desirability of each reaction are exactly equivalent to this matrix operation. Net influence of infinitesimal/incremental increase (NIII) -An accurate term for the OCARINA metric of node influence on dynamic VR graphs.
Here a more accurate term for the node-influence measure in [15] (when applied to a dynamic VR graph of the type shown in Figs. 2 and 3) is rationalised. The node-influence method in [15] returned the "demand" for a node based on the following logical propositions: If a relationship (R i ) positively affects X, this relationship is desirable. If a relationship (R i ) negatively affects X, this relationship is undesirable etc. However, the terms "desirability" or "demand" (a term used in [15], "electron energy demand"), though the latter is appropriate for some chemical processes (e.g. chemical oxygen demand), could be misleading when interpreting the output data of this node influence method on a VR graph. The algorithmic process is the sum of all the positive predecessor in-edges and all the negative predecessor in-edges of the target node. The data generated depends on how the edges are weighted, and in the case of the dynamic graph, where the derivatives are used to weight the edges, the result could be thought of as the sum of the positive and negative differentials, or "net differential". However, "differential" is already such a commonly used term that "net differential" could also be misleading when used in this context. Since a differential is the rate of change of one number with respect to another, which as mentioned above amounts to the influence of the independent variable on the dependent one, the term "net influence" would be more convenient. Given that such an 'influence' of a differential, by nature, is infinitesimally small, and that the change in the independent variable of a differential is always positive, and is thus an increase (from its current value to a higher value), and the influence itself is the sum of many individual influences, i.e. a net influence, an appropriate term might be "net influence of infinitesimal increase" or "relative net influence of infinitesimal increase" as the numbers are on a relative scale. A term like this should be of more assistance when interpreting results.
Clearly, making infinitesimal changes is not much use without an infinite amount of iterations, therefore the values of the net influences of infinitesimal increase that the algorithm produces, are best used to inform incremental changes as long as the following assumption is valid: In most cases this will clearly become less accurate the larger the difference is between the incremental increase and the infinitesimal increase. Thus, the smaller the increase the more likely it is to have the "desired" proportional effect on the target (as the gradient is less likely to change much over a short difference in value). In the same way, the gradient is also less likely to change much over short distances in the opposite direction, suggesting that decreasing the value of the variable if its influence is negative should have the same effect as increasing the value of the variable if its influence is positive. We can then amend the logical propositions in above to give actionable changes to be made to be made to the values of predecessors of the target node (Table 1) for future simulations, to increase the value of the target node. The following sections test this in an air plasma reaction set, on target nodes of O 3 , O and NO, focussing on electron energy as the key independent variable of influence.

Testing the Application of VR Graph Node-Influence Analysis Methods (Eigenvector Centrality, Katz Centrality, and OCARINA) in Time-Dependent 0D Plasma Chemical Kinetic Simulations
The VR graph and node influence analysis methods described in the previous section were investigated by applying them to a 0D plasma chemical kinetic simulation to analyse the net influence of the electron energy node on the formation of target chemicals. Sakiyama et al., [27] published a list of rate coefficients for a DBD plasma in moist air, with rate coefficients as functions of electron energy and temperature, and these rate coefficients were used for the 0D simulations in this work. Electron energy peaks of the type approximating Dielectric Barrier Discharge (DBD) streamers [28] were simulated by a triangular peaks increasing and decreasing linearly to and from a peak value. The gas temperature was kept It also follows that:

Net potential influence of incremental increase ≈ Net influence of infinitesimal increase
The greater the NIII of a variable (V i ) on X, the greater the change in X will result from a change in V i , and therefore the stronger the influence at a constant 300 K to limit the focus of the study to the effects of varying the electron energy. The simulation was run for 3,200,000 timesteps, with each timestep being 1 × 10 −13 seconds. Species concentrations and reaction rates were calculated for each time step and the results collected.
A dynamic variable-relationship graph ( Fig. 2) was constructed by adding relationships and variables from a list of reactions in [27] in csv format to a directional graph object in python using the NetworkX python library [29]. The differentials of the relationships were taken using the SymPy python library [30] and evaluated from values of the corresponding variables and relationships at any specified time step in the simulation. A section of an ego network of the electron energy (labelled ee) with a depth of 2 is shown in Fig. 6 and part of the full graph is shown in Fig. 7 using graph visualisations generated using Gephi [21]. Reverse graphs were also generated where the directions of all the edges in the graph were reversed. For an initial trial of the method, the O 3 variable node was used as a target node since efficient and/or high yield ozone formation in plasma discharges has been extensively studied for over a century. Thus the ideal electron energy condition sequence for ozone formation has already been established and could be readily compared with the results of the net-influence analysis.
The Eigenvector Centrality Index (ECI), Katz Centrality Index (KCI) and OCARINA/ Net-Influence of Infinitesimal Increase (NIII) were each generated for every 1000 th time step in a "default condition sequence" of a continuous series of successive 1 eV electron energy peaks. The ECI and KCI were both generated from the graph using functions eigen-vector_centrality_numpy() and katz_centrality_numpy() in the NetworkX library. Values of the graph edge weights needed to be scaled in order for the functions to work. A logarithmic scaling function was used: where n is the number to be scaled, and α is the scaling factor (in this case 1000). If n < 0 then: n scaled = 10 log 10 (n∕ ) 0 − n scaled = 10 log 10 (n∕ ) Fig. 7 The same ego network as in Fig. 6 but with the rest of the graph not hidden. The same size and colour criteria as Fig. 6 were used 1 3 was used. The OCARINA/NIII process did not require any scaling.
The results of this are shown in Figs. 9, 10 and 11, where any values are normalised, the normalisation function used was as follows: The results of the NIII of electron energy on O 3 formation were then compared to different electron energy condition sequences to assess how the NIII compared to the condition sequences that gave rise to greater ozone formation. Two different sets of sequences were compared. A set where the number of successive 1 eV pulses was varied from 1 to 10 (Fig. 8a) and a set where the relative peak values of two successive pulses were increased or decreased by 0.1 eV increments (Fig. 8b).

Node Influence Methods on a Default Condition Sequence of Continuous Electron Energy Peaks
The eigenvector centrality index did not appear to have any pattern (Fig. 9), this was likely due to the directional VR graphs containing nodes with no predecessors (in-degrees of 0). As such nodes will almost inevitably be present in VR graphs, the eigenvector centrality index is unlikely to be suitable for analysis of these graphs. The Katz centrality, however, solves this with the value beta, which is added to the centrality measures from the very start of the calculations, essentially giving each node a very small amount of importance to prevent nodes from being assigned a centrality of 0 even when they have high out-degree measures. For this reason there is a much clearer pattern in the Katz centrality results in Fig. 10. The Katz centrality of the electron energy node did not seem to vary much over the course of the simulation, however, the reverse Katz centrality (the Katz centrality on the reversed graph) of electron energy appeared to decrease sharply over the course of the first simulated electron energy pulse, and then again at the start of the second. This could be explained by the rapid rates of electron impact reactions resulting from the increase in electron energy during the first peak, giving rise to increased reaction rates,  and hence edge weights. These initial electron impact reactions then produce products that undergo many further reactions that are independent of electron energy, hence the relative importance of the electron energy decreases. This suggests that the reverse Katz centrality index is representative of the influence of electron energy on the whole graph, and may be an appropriate method to use for assessing the change of node influence with time for VR graphs.
The NIII was evaluated at three preceding depths of variable nodes from a target node of O 3 . These NIII values appeared in either positive or negative peaks which coincided with the peaks of electron energy (see Fig. 11). Interpreted using Table 1, this suggested the following: For electron energy as a depth 2 ("direct") predecessor to ozone, a positive change in the peak value of the first electron energy pulse would have a small net negative change in the value of the ozone concentration. After this first pulse, positive changes in the peak values of the electron energy pulses would have large net negative changes in the value of the ozone concentration.
For electron energy as a depth 4 predecessor to ozone, a positive change in the peak value of the first electron energy pulse would potentially have a relatively large net positive change in the value of the ozone concentration. Any positive change in the peak values of subsequent electron energy pulses would have increasingly less positive change in the value of ozone concentration, and this eventually levels off. Conversely, any equal and opposite changes to any of these pulses would have the equal and opposite effects. Results of the OCARINA node-influence ranking algorithm of electron energy on influence on ozone formation over time in a 0D plasma chemical kinetic simulation of an air plasma, compared with the electron energy over time. The simulation simulated 10 ns triangular pulses of electron energy to simulate the effects of plasma streamers. Three depths of predecessor were analysed, direct precursors of O 3 (depth 2 predecessors) precursor of direct precursors of O 3 (depth 4 predecessors) and precursors of precursors of direct precursors of O 3 (depth 6 predecessors). Due to the differences in the orders of magnitude of the different precursor phases, the results were normalised to allow comparison of the trends I.e. decreasing the peak values of each successive electron energy pulse would have a less strongly negative effect. The trend in the NIII at depth 6, and its implications are the same as that of depth 4 in this case. Clearly, whilst depths 4 and 6 agree with each other, they are contradictory to depth 2, in that they suggest that increasing the electron energy is overall positive to O 3 formation, whereas depth 2 suggests the opposite. This discrepancy was fortunately resolved in the results interpretation of simulations with the condition sequences shown in Fig. 8b, which will be discussed in the next section.
There is also some agreement between the NIII at these three depths, however, in that the first electron energy peak is the least negative/most positive to O 3 formation. This is also reflected, at least to some extent, in the reverse Katz centrality which suggests that the initial electron energy peak has the most influence on the whole network, and successive peaks less so. In contrast to the sharp decreases of the Katz centrality of electron energy of the first two peaks, the NIII of electron energy on O 3 at depths 4 and 6 appear to show a more gradual decay in the positive influence of increasing the electron energy on O 3 formation. This suggests that each successive electron energy pulse may be giving diminishing returns for increasing the ozone concentration.
To test if this was indeed the case simulations were run with increasing numbers of electron energy pulses from 1 to 10 (see Fig. 8a), and each simulation was then left to run for a total of 3.2 × 10 −7 s. The final O 3 concentrations of each simulation are plotted in Fig. 12, along with the final O 3 concentration per electron energy pulse as a measure of energy efficiency.
It can be seen in Fig. 12 that each additional electron energy pulse does indeed result in a progressively smaller increase in the ozone concentration in a trend strongly correlating with the first 10 "peaks" of NIII depths 4 and 6 of electron energy on O 3 in Fig. 11. The energy efficiency measure in Fig. 12 suggests a single pulse would be 336% as energy efficient as 10 successive identical pulses. This suggests that correct interpretation of the NIII results of a single "baseline" simulation, could anticipate a superior conditional sequence of pulses of 1 eV, 0 eV, 0 eV, 0 eV, 0 eV, 0 eV, 0 eV,… etc. for maximum output per input, Fig. 12 Results of the final ozone concentrations and final ozone energy efficiencies of formation in ten different simulations of incrementally increasing numbers of streamers (from one to ten, see Fig. 8.a) 1 3 and a superior conditional sequence of 1 eV, 1 eV, 1 eV, 1 eV, 1 eV, 1 eV, 1 eV,… etc. for maximum output. The results also suggest that increasing the peak values of any of the pulses beyond 1 eV would also increase the O 3 formation, and decreasing it would have the opposite effect. However, Fig. 12 does not show the effects of all possible condition sequences of 10 pulses of peaks of either 1 eV and 0 eV, or of peaks of any other values, so it cannot be concluded that the NIII is indicative of the ideal condition sequence. Testing this is much more challenging. Testing all permutations of 10 pulses of either 1 or 0 eV, would require something on the order of 2 10 simulations, and increasing and decreasing the each pulse by only a single increment either way 3 10 , much less any more peaks/ increments.

Comparison of NIII for O 3 , O and NO Targeting with Results of Incrementally Increasing or Decreasing Two Successive Peak Values by 0.1 eV
For an initial look at this, a preliminary study of just 2 pulses increased and decreased by a single increment gave a much more manageable number of just 3 2 different condition sequences to simulate (Fig. 8b). Simulations of these 9 different condition sequences were run for three different chemicals, O 3 , O and NO. The gas temperature was also varied by ± 20 K coincident with the same time intervals as the electron energy pulses, but these increments were too small to make any significant difference in any of the outcomes, and constraints on time did not permit larger increments to be tried. For this reason, the NIII of gas temperature on the concentrations of O 3 , O and NO were not included in this work, and a constant gas temperature of 300 K was used throughout. The NIII's for each target chemical were generated from a "baseline"/default condition sequence of two consecutive 1 eV pulses.
The results of these simulations are shown in Figs. 13, 14, 15 and 16, input (electron energy in this case), productivity (value of targeted output, species concentration of either NO, O and O 3 in this case) and cumulative efficiency (value of targeted output/cumulative input) are shown for each targeted chemical and compared with NIII of electron energy for each species for three depths of predecessor. There was good agreement between the NIII results and the results of the condition sequences for each of the target chemicals. Additionally, the different NIII depths were found to correspond to different time phases in the simulations.
The NIII of electron energy (ε) on ozone (Fig. 13), for example, is (a) negative for the direct precursors (precursor 1/predecessor depth 2) for both electron energy peaks, (b) more strongly negative for the second peak than the first. This suggests that (i) the effect of an increase in either of the two electron energy peaks would result in a decrease of ozone formation (and vice versa for a decrease) and (ii) according to the NIII for depth 2, an increase in the second peak would decrease ozone formation more strongly than an increase in the first. Point (i) is supported given that there is no O 3 formation during the rise of the first peak, and the ozone formation then increases as the first peak falls (Fig. 13). Ozone formation over the course of the second peak is likely to be due to the effects suggested by the NIII's of depths 4 and 6, as there will inevitably be an overlap in the effects of predecessor depths 2, 4, and 6 once O 3 begins to be formed. It seems that as soon as the instantaneous electron energy ceases (at 0.2 × 10 -7 s and beyond) the effects "predicted" by depth 4 and 6 predecessors can be immediately observed, which suggests that the effects of NIII's of depth 2 correspond to an instantaneous influence on the target node (in a network of this type).

3
The depths of 4 and 6 NIII's of electron energy on ozone were both (a) positive for both peaks and (b) more positive for the first peak than the second peak. This suggested that (i) increasing either peak should increase the ozone formation and (ii) increasing the   first peak should have a stronger effect than increasing the second. Support for point (i) can be observed from the fact that the lowest decrease of both peaks (0.9 eV,0.9 eV) gave the lowest ozone concentrations, and the highest increase of both peaks (1.1 eV,1.1 eV) gave the highest (Fig. 13). Point (ii) can be supported by the observation that, for the same total energy input, a higher peak followed by a lower peak, invariably corresponded to a higher ozone formation than a lower peak followed by a higher peak in every case (compare 1.1 eV,0.9 eV vs 0.9 eV,1.1 eV for example). The fact that both predecessors 4 and 6 are positive corresponds with the continuing increase of O 3 , although since predecessor depths 4 and 6 are both positive and have almost identical relative magnitudes in their peaks, it is not possible (at least within the scope of this work) to tell where the effect of one begins and the other ends, and the nature of the overlap between them, though the two effects would likely compound.
Some similar agreement between NIII and condition sequences can be observed for O targeting (Fig. 14). The depth 2 predecessor values are generally positive for NIII of ε on O which seem to correspond to a very strong instantaneous increase in the O concentration. The second increase being a smaller increase than the first, also corresponds to the second NIII "peak" being smaller than the first.
After the initial electron energy pulses, the subsequent trends in O concentration correspond with the NIII results for predecessors 4 and 6. The results of all simulations appear to increase up to approximately 0.8 × 10 -7 s, corresponding to a positive depth 4 predecessor for both peaks, and a subsequent decrease for all simulations after ~ 0.8 × 10 -7 s corresponds to a negative depth 6 predecessor for both pulses. The ratio of the NIII for the two electron peaks also corresponded with observations from the different simulations, as in a similar manner to the O 3 results, the electron energy pulse has a higher NIII than the second, and the same trend (higher first peak giving a higher target chemical concentration than a higher second peak for the same total energy input) was observed. The fact that depth 6 predecessor results had a more strongly negative first peak than the second peak, suggests that running the simulation for long enough might see these trends reverse. It would indeed be interesting to see if continuing the simulation for longer would result in the O concentration of the 1.1 eV, 1.1 eV simulation eventually decreasing below the 0.9 eV, 0.9 eV simulation, and 1.1 eV, 0.9 eV decreasing below the 0.9 eV, 1.1 eV, but unfortunately this was outside the scope of this study.
The NIII's for NO (Fig. 15) appeared to have two broad peaks that coincided with the two electron energy peaks which appeared to have similar relative magnitudes to the O and O 3 results, and also two very narrow peaks coinciding exactly with the ends of the electron energy peaks. For NIII depth 2, the two narrow peaks had very similar magnitudes, with the first peak having a slightly greater absolute value. The two narrow peaks were absent in NIII depth 4 and present in NIII depth 4 and both were negative with the first peak being between 4 and 5 times the absolute value of the second. The very sharp drops seen in the depth 2 predecessor results of the NIII also appear to correspond with a sharp transition in the concentrations at these points in time.
The broader peaks for NIII depth 4 and 6 appear to show agreement with the simulations in a similar way to O and O 3 : increase following the initial electron energy peaks (corresponding with positive depth 4 NIII), then a slight decrease (corresponding with negative depth 6 NIII). There is also the same agreement as O and O 3 with the different combinations of peak heights in the simulations, and the NIII results (higher overall resulting in more NO formation, and higher then lower peaks forming more NO than lower then higher peaks for the same total electron energy). The very steep straight line dips could possibly be due to at least one reaction in the rate coefficients (reaction e + O 2 O 2 − [27]) having a step function relationship with two different rate coefficient functions for different electron energy levels. The methods used in this paper currently do not account for this, and future modifications will need to take this possibility into account.

Discussion
Everything in the previous sections suggests a relationship between the NIII and the variation in output resulting from changing the inputs in the simulation over the ranges investigated. Furthermore, the results for O 3 , O and NO (Figs. 13, 14 and 15) all suggest that, for the most part, each antecedent graph-theoretical depth from the target node appears to correspond with each subsequent temporal phase in the simulation, resolving the issue of the contradictory NIII results for the different depths mentioned in previous sections. This, and the correlation of the results in Figs. 11 and 12 suggest that the OCARINA method can be effective for targeted condition sequencing for these 0D time-dependent study. More broadly, the results suggest that more locally optimal conditions for targeted chemical outcomes can be identified by the OCARINA method of analysis of the relationship between variables in the context of a network. This has been achieved by a variable-relationship bipartite graph such as the one described in Figs. 1, 2, 3, 4, 5 which, once made, can be evaluated at any point in a time-dependent simulation by a simple matrix operation on the graph adjacency matrix. The OCARINA method should even be scalable to multiple spatial elements, scaling the method up to 1D and higher, provided this can be represented by a variable-relationship graph. With finite element methods, for example, once a finite element array is defined the numerical approximation to the PDE (Partial Differential Equation) system is a network of nodal ODEs (Ordinary Differential Equations) [31]. The variable relationship graphs for this would be much larger than those used in this work, as an increase in the number of spatial elements would mean an exponential increase in the number of elements in the adjacency matrix of the graph, and an increase in the number of dimensions often means an exponential increase in the number of spatial elements. Fortunately, the simulations would also scale in a similar manner, thus the OCARINA method should scale in the same way as the model it is used on. Given that the main components of this analysis are matrix operations, the OCARINA method could also be shared across multiple processors if the application needs to be scaled up for more complex systems.
Since the Katz centrality also operates on the adjacency matrix, this would also scale accordingly. The analysis of the results in Fig. 10 suggest that the Katz centrality gives a measure of the relative influence of a node (in this case electron energy) over time with respect to the whole graph. Although the Katz centrality is described as a solution to eigenvector centrality for directed graphs, there were still issues with its use. As discussed in [32], there is not a strict definition of what makes a centrality measure, and researchers are constantly proposing variations to existing centrality measures. However, Lawyer [33] suggests that while centrality measures are good at highlighting the very most influential nodes, they do not do a good job of quantifying the importance of the other nodes (midimportance nodes for lack of a better phrase). Identifying the most influential node in a network is useful, as in the absence of any other information, this might be the node that is most worth attention. Nevertheless, it is very unlikely that this node will have equal influence on all the other nodes, and may have little influence on a particular target node. Thus Katz centrality values on their own are not sufficient for targeted condition sequencing.

Limitations of the Methods in this Study
Whilst the OCARINA should be scalable to be applied to more complex models and higher dimensions, the difficulties of verifying the OCARINA algorithm results also scale accordingly. To be as rigorous as possible would require all possible condition sequences to be investigated and compared to the output of the OCARINA algorithm. Ideally a parametric sweep of a range of different values of all ten pulses would have been simulated, but even as little as ten values of each pulse (0.1 eV, 0.2 eV, 0.3 eV … 1 eV, for example) would result in 10 10 combinations, each one needing a separate simulation. Unfortunately, this was computationally unfeasible within the timescales of the project and the available computing power. For this same reason, NIII's of multiple variables were not looked into (except for the omitted gas temperature results mentioned above).
Another limitation of the study was that the main application of the algorithm was on the electron energy. Whilst this is clearly the right variable to focus on for plasma chemical waveform targeting, it is also a variable with no predecessors. This study is therefore limited to observing the effects of altering completely independent variables, and not variables which are embedded in the cycles of the network (such as most of the reactions and species) and thus are both influenced by as well as influencing the target node. This also includes the target node itself.
Some additional issues also need to be discussed. Looking at Figs. 1, 2, 3, 4, 5, 6, and 7, one might wonder how increasing one of the variables might inadvertently affect the others. Clearly the weight of each variable to reaction edge depends on more than just the start and end nodes, they are functions of the other variables in the relationship. In Fig. 3, for example, changing the electron energy would change other things outside of this reaction that would indirectly change e and O2 concentrations, and this would alter the effect that increasing the electron energy would have. This might mean the very means of defining the effect of increasing the electron energy (NIII), change as soon as it's increased. Whilst this is true, this is why the term "Net Influence of Infinitesimal Increase is accurate, as it is only strictly true if every variable is changed infinitesimally. As mentioned previously, applying the OCARINA method to the dynamic graph works under the assumption that an increase in small enough increments will be a close approximation to the infinitesimal increase. Also, as shown in Figs. 13, 14 and 15, the effects of a node on all other nodes and on itself is effectively increasingly taken into account as the node depths of analysis increase. Each depth of analysis is done on the whole adjacency matrix, including the target node. The results in Figs. 13, 14 and 15 also suggest that indirect "feedback" effects of a node on itself are not instantaneous, as this is a time dependent simulation, only a certain amount of change can happen in a given time step. This being the case, it would be interesting to include 'time' itself as a node in these variable relationship graphs in future. Figure 12 shows that it is possible in some cases to make inferences outside the small ranges of the small variable space covered by the 9 different conditional sequences in Fig. 8b, and whilst results like this cannot be guaranteed, it suggests there are scenarios in which changes can be made in greater increments than those shown in Figs. 8b and 13, 14, 15 and 16. It would be convenient if these scenarios could be identified, and it may be possible to do this using VR range graphs, but this will have to be investigated in future studies.

Estimating the Increments
Being as precise as possible, the y-axis NIII plots ought to be labelled "Predicted net relative change in the value of the target node as a result of an infinitesimally small positive change in the value of the given source node at its current value". For a dynamic graph, all these numbers are true only at the immediate values of the source nodes. Knowing how much to change the inputs by is key to being able to make further use of this novel method. If the NIII results indicate that increasing an input would be strongly conducive to increasing an output, should the input be increased by 10%? 20%? 50%? Etc.? What if increasing it gives the opposite result?
The results of the NIII are true for an infinitesimally small increase, essentially a single point. The closer a real increase is to this single point, the closer the result of the increase will be to what the NIII suggests it will be. Thus the smaller the increase the better, and an algorithm that moves over very small steps will likely be the most accurate. This will of course be more time-consuming, so increasing the step sizes in changes in the variable will be desirable. However, too great an increase may overshoot the peak in the relationship and begin a decrease, and so identifying where the peaks are is the missing ingredient to enhancing this method. Further to the previous point, the NIII is composed of a large number of derivatives, which are all only true for single points along the graphs of each relationship, and the distance from each point to a local maximum (peak) of each relationship must be different. It is possible to know where each differential is on the x and y axis in each relationship, and therefore how far away from the local peak that is. Local maxima and minima of each relationship can be found before any simulations are run and stored as edge attributes, and the difference between those and the point of the NIII for each relationship could be measured at any point during a simulation. This could be done for each edge and the resulting values used for an error measurement for the NIII. Local maxima and minima will vary depending on the values of all the independent variables for each relationship, so 2, 3 or higher dimensional shapes might be needed to represent this error (which could be represented as an n-dimensional identity matrix, for example), which again could be calculated before the simulation is run. This could point the way to assessing how large an incremental increase or decrease should be. Finally, it must be noted on making incremental changes that if a variable is changed earlier on in the simulation, this will effectively invalidate the NIII readings after that point. It may be best to make changes to the last time steps first.

Potential Variations of Set 0: Multiple Targets, Variable Targets and Anti-Targets
Looking at Set 0 The only value of this set used in the example was D 0 = 1, all the other variable values = 0. Clearly, if D 0 was set to -1, and the same method run, then the output of the algorithm would be conducive to reducing the value of D 0 . Hence D 0 would be an "anti-target". This being the case, the algorithm could be run with one target, and one anti-target. Say, for example, we want to avoid a troublesome by-product of a process, we could set this as an anti-target. If there are any pathways that favour the target without the anti-target, the NIII will be composed of these, as any pathways shared by the target and the anti-target will be cancelled out.
It also follows that more than one variable can be set as the target (if two chemicals are required for example) and that they could be given different set 0 values if there is an inclination to produce them in different ratios. Clearly there will be many scenarios where this will not be possible due to the number of pathways shared by the multiple targets/antitargets. A preliminary graph analysis could be run to assess the number of shared relationships/pathways to inform whether such a multi-target analysis would be worth pursuing.

How this Translates to Systems in Reality
Using these graph theory methods are only as good as the models that the graphs are formed from. Anything that applies to the model, applies to the graph. If the models do not accurately represent reality, then the graphs will not either. If big assumptions are made to simplify the model, these will apply to the graph as well. Assuming the model is an accurate representation of a system in reality, the results of this algorithm highlight changes which can be made to the variables, which correspond to actions in the real system.
These results are based on simple electron energy pulses aimed to roughly simulate streamers in a DBD plasma source. The relationship between the voltage waveform for this plasma source and the DBD streamers is not a direct one and is dependent on the frequency and the plasma source itself. The methods described in this paper, however, could be applied to a mathematical model that accurately represents the relationships between the power supply, voltage waveform, plasma source and DBD streamers, assuming such a model is practical. In the absence of such a model, the NIII of the electron energy could be used to identify the ideal sequence of electron energy magnitudes, which in practice could be supplied either successively in time, or successively in space (e.g. if a gas flow is moving through successive discharge conditions). Aside from electron energy and gas temperature, if a reaction relationship, for example, has a high NIII, an action would be to somehow catalyse that reaction. If a chemical variable, for example, has a high NIII, an action could be to inject some of that chemical into the system at the time the NIII is highest.

Wider Applicability
More generally the question of how much electron energy to supply is a question of when to supply any input for a given output. Put in the simplest terms, the question of "how much of what to supply and when?" is an almost universal one. The answer to this is "when it matters most", and something could be said to matter most at the point it has the strongest influence on an intended outcome, i.e. where it will make the biggest difference with a capital Δ.
For any system that can be adequately described by a series of time-varying equations, a variable-relationship graph can be made by forming edges and nodes from each of its equations as above to form a weighted directional bipartite graph containing both positive and negative edges. As a great number of systems can be described by such equations, the techniques outlined in this work could have much wider applicability not just in reaction engineering but even in engineering in general. A general application of this work may therefore be a useful addition to the toolkit of techniques for improving the processes we all depend on.

Conclusions
It has been demonstrated that a plasma chemical reaction system, including electron energy and gas temperature, can be represented as a dynamic bipartite variable-relationship (VR) graph. The results suggest that the OCARINA method of analysis can be used to identify more effective conditional sequencing of electron energy for a target chemical species.
Using this dynamic VR graph, the OCARINA method described in [15] can be used at any point on a time-dependent 0D chemical kinetic simulation, and can also be simplified to a vector-matrix product of the VR graph adjacency matrix with a target identity vector. This matrix operation can be repeated for successive 'depths' of predecessors of the target node without any loss of computational expediency. The resulting set of values give a measure of the 'net influence' of increasing the value of each node on a target node, which are different at different depths.
0D chemical kinetic simulations were run with 1-10 electron energy pulses and it was observed that each successive pulse gave diminishing returns on ozone output. This agreed with the 'net influence of incremental increase' (NIII) results of running the OCARINA algorithm for a baseline of continuous pulses. The eigenvector centrality index (ECI) and Katz centrality index (KCI) were also run the graph for this same baseline simulation and it was found that the ECI was not suitable for this type of graph, but the reverse KCI gave meaningful values. 0D chemical kinetic simulations were also run for each combination of ± 10% of the peak values of just two electron energy pulses. OCARINA was run on a baseline of the default values of these two pulses. The results of these kinetic simulations largely agreed with the NIII results of the OCARINA algorithm for target nodes of O 3 , NO and O. It was also observed that the NIII for each antecedent depth of predecessor, appeared to correspond with a successive temporal 'phase' in each simulation.
Simple alterations to the target identity vector could also allow OCARINA to be run on other target species, and even on multiple target chemical species and also negative target nodes, though this will need to be examined in future work. These methods could also potentially scale to 1D, 2D and higher dimensional models, and could potentially be applied to any time dependent set of equations, provided they can be represented with a VR graph.
Future work should therefore include investigating the effects of altering variables other than the electron energy, in finding the limits of how these variables can be increased in response to the NIII results (i.e. to more rapidly approach optimal conditions) and also in applying these methods to other systems and higher dimensional models. Further development of these methods could be a key component of a dynamical systems alternative to the steady state, mixed-integer nonlinear optimisation method.

Appendix
Here the steps are outlined for the generalisation of the logical statements in the OCA-RINA method [15] to matrix operations on variable-relationship graphs (Figs. 1, 2, 3 and 4) to verify that the two are essentially equivalent, and that the matrix operations can therefore be applied as a node-influence ranking method.
To begin, a target node (X) must be specified in the variable relationship graph (in this case a variable node) which represents the outcome we wish to change (if we are targeting ozone for example, the target node would be the O 3 variable node). Given that any node in a directional graph can only be directly influenced by a node directly preceding it, we need only be concerned with predecessors. Given also that a variable is always preceded by a single relationship, which in turn is always preceded by one or more variables (unless the variable is a completely independent variable), any node affects the target node through a series of successor nodes as follows: The total of all such pathways between any node and X constitute the total/net influence of this node on the target node X.
Clearly a node could be present at multiple points (depths) along a single chain and be present at different depths along different chains. It therefore follows that the node influence will be different at different depths away from the target node (X). This being the case it is important to be able to evaluate the net-influence of each node at different depths. This can be done using the logic in Tables 2, 3 and 4, and using the relationships in the following sections.
The influence of a direct predecessor relationship (R i ) on the target node X can be assessed using the logic in Table 2. In turn, the influence of a variable V 1 on X can be assessed by assessing its influence on R i , given R i 's influence on X, if V 1 is a direct predecessor of R i (Table 3). Going a further step back through the predecessors in the graph, this logic can then be applied again as shown in Table 4, and so on. In such a graph, the simplest measure of the "net influence" of any node on the target node, is the sum of all of its pathways of influence on the target node X (all positive pathways minus all negative pathways) at a given depth.
Nodes must be grouped into one or more sets depending on if they satisfy one or more of these criteria. Set O is the set containing only the target node. All other sets are sets of predecessors to the target node. Each set is named according to how its elements affect their successors. Set PN, for example, contains only nodes which positively affect nodes which negatively affect the target node, whereas set PNPP contains only nodes which positively affect nodes which negatively affect nodes which positively affect nodes which positively affect the target node.
Each node in each set is assigned a value to measure the extent of its effect on the node X. The same node can be present in more than one set, but its value calculation will be specific to that set. These values of each node (species or reaction) in each set are evaluated using the simplest operations possible, multiplication and addition, as follows:  is desirable If a relationship (R i ) positively affects X, and a variable (V 1 ) negatively affects this relationship, this variable (V 1 ) is undesirable If a relationship (R i ) negatively affects X, and a variable (V 1 ) positively affects this relationship, this variable (V 1 ) is undesirable If a relationship (R i ) negatively affects X, and a variable (V 1 ) negatively affects this relationship, this variable (V 1 ) is desirable If a relationship (R i ) positively affects X, and a variable (V 1 ) positively affects this relationship, and a preced- If a relationship (R i ) positively affects X, and a variable (V 1 ) negatively affects this relationship, and a preced- If a relationship (R i ) positively affects X, and a variable (V 1 ) negatively affects this relationship, and a preced- If a relationship (R i ) negatively affects X, and a variable (V 1 ) positively affects this relationship, and a preced- If a relationship (R i ) negatively affects X, and a variable (V 1 ) positively affects this relationship, and a preced- If a relationship (R i ) negatively affects X, and a variable (V 1 ) negatively affects this relationship, and a preced- If a relationship (R i ) negatively affects X, and a variable (V 1 ) negatively affects this relationship, and a preced- As every value in every set is effectively a multiple of the value given to X in set 0, the default value of X in set 0 is 1.
The "desirability" or "demand" of a node (for increasing the value of the target node) is the sum of all the parallel values of a node at the specified number of predecessor steps from the target node. A variable node (V) for example, would have its "desirability" ranked/scored as follows:  16 A map of the logic of the predecessor ranking, the predecessors are ranked according to how they influence a specified target node. This figure also serves to illustrate the problem of computing this logically, that as more precursive stages are introduced, the number of sets, and therefore the number of computations, effectively doubles Nodes with net positive or negative ranks should result in an overall increase or decrease in the value of the target node respectively.
The dependencies of these calculations on the sets are delineated in Fig. 16. It is clear from Fig. 16 that with each additional depth of predecessors to the sequence, the number of sets that would need to be processed doubles. This would lead to predecessor ranking becoming computationally unfeasible this way if the number of predecessors is too great.
This raises the question: If the process was done by adding the sets together at each depth before evaluating the next depth, instead of evaluating a negative and positive set for each set in the succeeding depth, would this produce the same results? Or would some pathways of influence cancel each other out if added together, which would eliminate effects at preceding depths being taken into account?
Conveniently, it can be shown that adding the sets together at any depth to make a single set, and evaluating preceding sets starting from that, is exactly equivalent to the above process, and thus won't change the outcome at any depth. The following section verifies this.
Reduction of the set summation operations to adjacency matrix operations. Using the graph in Fig. 4 as an example, since D is the target node, we set the value of D as 1, and all the other values as 0 in set 0.
Target node set: Each node "score" in set P and set N are then functions of the values in set 0, as follows: Depth 1 predecessors: Each node "score" in set PP, NP, PN and NN are then functions of the values in sets P and N as follows: Depth 2 predecessors: And for all sets in the next stage: Depth 3 predecessors: Relationships: And for all sets in the next stage: Depth 4 predecessors: Variables: Expanding the brackets means: Depth 4 predecessors: T e4 = T ePPPP + T eNPPP + T ePNPP + T eNNPP + T ePPNP + T eNPNP + T ePNNP + T eNNNP + T ePPPN + T eNPPN + T ePNPN + T eNNPN Depth 3 predecessors:

Depth 2 predecessors:
The relationships in the depth 1 predecessor stage (succeeding the variables in the depth 2 predecessor stage) are either positive or negative, therefore either A PP + A NP or A PN + A NN will be zero, and hence the value of the variables in the depth 2 predecessor stage will be: Depth 2 predecessors:

Depth 1 predecessors:
Calculating all phases then simplifies to: Target array: And the Final rank set would then be defined as follows: Depth 4 predecessors: A 4 = r 3 R 13 + 0 R 23 + r 6 R 33 This is equivalent to a series of matrix operations, which can be performed on the adjacency matrix of the graph, the values of all nodes in each set being calculated as follows: Adjacency matrix: To get depth 1 predecessors: To get depth 2 predecessors:

To get depth 3 predecessors:
To get depth 4 predecessors: The above matrix operations can be represented in the following notation: where x i is the ith column of the node ranking vector, a i,j is the jth row of the ith column, and n is the depth of predecessor from the target node. The initial node ranking vector (n = 0) with the index of each node set to 0, except the target node which is set to 1.
All the above is effectively a simplified generalisation of the OCARINA (Optimal Condition Approaching via Reaction-In-Network Analysis) algorithm in [15] using matrix operations. Taking account of the generalisation into a Variable-Relationship graph, the same acronym can still be used, with a slight modification: Optimal Condition Approaching via Relationship-In-Network Analysis.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.