Problem formulation

where c _ij denotes the coefficient associated with the directed edge V _j  → V _i , and ε _i denotes the disturbance error term that follows a certain distribution (Gaussian or non-Gaussian [40, 41]) with mean zero. For simplicity, all disturbance error terms are assumed to be independent. By definition, E specifies the structure of the coefficient matrix C = [c _ij ], i.e., c _ij  = 0 if no edge exists in E from V _j to V _i for i ≠ j. When a network structure contains one or more loops, G is a directed cyclic graph (DCG) and the corresponding SEM is called a non-recursive model; otherwise, G is a directed acyclic graph [42] and the corresponding SEM is called recursive. Although Drton’s condition [31] and Foygel’s half-trek criterion [32] are applicable to the identifiability analysis of non-recursive SEMs, the identifiability of parameters on a loop may be still inconclusive. Due to the lack of mature structural identifiability analysis techniques for examining every unknown parameter of a non-recursive SEM, this study focuses on recursive SEMs (i.e., DAGs) only.

Definition 1 (observation strategy). Given a graph G = (V, E), its observation strategy can be denoted by a binary vector \( O={\left({O}_{V_1},\cdots, {O}_{V_n}\right)}^T \), where \( {O}_{V_i}=1 \) if node V _i is observed and \( {O}_{V_i}=0 \) if V _i is unobserved.

where \( {\displaystyle \sum_{i=1}^n{O}_{V_{{}_i}}} \) is the total number of observed nodes in an observation strategy O, and 1 denotes a vector of ones. For clarification, we stress that the observation measurements are for the random variables associated with network nodes, and we assume (n – m) of them are observed in the original observation strategy, where n denotes the total number of nodes and 0 < m ≤ n.

The objective function above is minimized with respect to the originally unobserved nodes, subject to the constraint D = 1. During the minimization process, it needs to be repeatedly verified whether all parameters have become identifiable (i.e., D = 1). For this purpose, an efficient algorithm for structural identifiability analysis of SEMs is needed. Here we consider the identifiability matrix method proposed by Wang et al. [34]. Briefly, structural identifiability of parameters can be verified by examining the number of solutions to the symbolic polynomial identifiability equations generated by Wright’s path coefficient method [43, 44]. To avoid the expensive symbolic computation involved in reducing such identifiability equations, the identifiability matrix method proposes to derive binary matrices from symbolic polynomials and thus enable us to determine the number of solutions via several simple matrix operations. It is noteworthy that Wang’s work [34] does not explicitly handle colliders involving bidirectional arcs when generating identifiability equations with Wright’s method, however, the identifiability matrix method is still applicable here as we do not consider bidirectional arcs in DAGs.

Identifiability gain and must-be-observed nodes

The optimization problem in the previous section is combinatorial in nature. Therefore, if the number of the originally unobserved nodes (denoted by m) is not small, enumerating all the 2 ^m different possible observation strategies over these nodes will be computationally expensive. We thus need an efficient algorithm such as dynamic programming to obtain the solutions. For this purpose, a few more definitions need to be introduced first.

Definition 2 (redundant identifiability equation). Given a set of identifiability equations, an identifiability equation IE(V _i ,V _j ) is redundant with respect to that set if it can be expressed as a linear combination of the equations in that set.

Definition 3 (cardinality of observation strategy). Given an observation strategy O for a network G, one symbolic polynomial identifiability equation can be generated for each pair of d-connected [28] observed nodes using, e.g., Wright’s path coefficient method. Then the total number of non-redundant identifiability equations is called the cardinality of O, denoted by f(O).

where the calculation of f(O) is a key challenge because it depends on specific network structure and observation strategy and thus has no closed-form solution. We thus introduce the following definition.

Definition 4 (identifiability gain). Given a network G = (V,E), let O ^(k) and f(O ^(k)) denote an observation strategy and its cardinality, respectively. Let V _i be an unobserved node in O ^(k), and only V _i becomes observed in a new observation strategy O ^(k + 1) with the observation statuses of other nodes remaining unchanged. Let f(O ^(k + 1)) denote the cardinality of O ^(k + 1). Then the identifiability gain of observing V _i , denoted by g(V _i , O ^(k)), is calculated as g(V _i , O ^(k)) = f(O ^(k + 1)) − f(O ^(k)).

By definition, g(V _i , O ^(k)) is the difference in cardinality between two consecutive observation strategies O ^(k) and O ^(k + 1). That is, after V _i becomes observed in O ^(k + 1), we need to find out the number of newly added non-redundant identifiability equations. First, if another node V _j (i ≠ j) is observed in both O ^(k) and O ^(k + 1) and there exists a Wright’s path [46] of length 1 connecting V _i and V _j , it can be shown that the newly added identifiability equation, denoted by IE(V _i ,V _j ), is non-redundant (see Lemma 1 and Additional file 1 for theoretical justification). However, if the length of every Wright’s path between V _i and V _j is greater than 1, the identifiability equation IE(V _i ,V _j ) is not always redundant, and it depends on both the node’s observation status and the structure of the network. Here we introduce the concept of detour-path before we further elucidate the redundancy issue. Consider a DAG G = (V,E) and two d-connected observed nodes V _i and V _j . Assume that there exists a Wright’s path P _ji between V _i and V _j as well as an observed node V _k (k ≠ i, j) on P _ji , and the direction of P _ji is from V _i to V _k and then to V _j . Now let P _ki and P _jk denote the two segments of P _ji , then P _ki entering node V _k has an arrow pointing into V _k while P _jk exiting node V _k has an arrow pointing away from V _k . However, if there exists another Wright’s path between V _k and V _j , denoted by \( {\tilde{P}}_{kj} \), which has no any other observed nodes besides V _k and V _j and has an arrow pointing into V _k , then V _k is a collider with respect to P _ki and \( {\tilde{P}}_{kj} \). Thus, we call the Wright’s path segment P _jk the detour-path, and call V _i , V _j and V _k the upstream node, the downstream node, and the collider node of the detour-path P _jk , respectively. By definition, a detour-path can have only one downstream node and one collider node but may have one or more upstream nodes. Moreover, multiple detour-paths can share the same upstream node, the same downstream node or the same collider node. Several examples are given in Fig. 1 to illustrate the concept of detour-path.

In addition, when an upstream node V _i is shared by two or more detour-paths that have the same downstream node, V _i is called a shared upstream node; otherwise, V _i is called an exclusive upstream node. Note that a detour-path can have both exclusive and shared upstream nodes in the same time, and the collider node of one detour-path can be an upstream node of another detour-path. Consider two detour-paths that have no exclusive upstream nodes, if they share the same downstream node and at least one upstream node, or one upstream node of one detour path is the collider node of the other detour-path, then two detour-paths are intersecting. One can tell that if \( {P}_{j{ k}_1} \) intersects with \( {P}_{j{ k}_2} \) and \( {P}_{j{ k}_2} \) intersects with \( {P}_{j{ k}_3} \), then \( {P}_{j{ k}_1} \) also intersects with \( {P}_{j{ k}_3} \). Then we consider a downstream node V _j , let S_IDP denote all the intersecting detour-paths, and let S_SUN denote all the shared upstream nodes of S_IDP. Similar to a single unknown parameter, the coefficient product \( W P={\displaystyle \prod_{e dg{ e}_l}}{\theta}_l \) of a Wright’s path P can be deemed as a single parameter and one can tell its structural identifiability based on identifiability equations. If a detour-path P has at least one exclusive upstream node, then the Wright’s coefficient WP of P is globally identifiable (see Lemma 2 and Additional file 1 for theoretical justification). Also, for a group of intersecting detour-paths, if the node number of S_SUN is equal to or greater than the number of intersecting detour-paths in S_IDP, then the Wright’s coefficient of each detour-path in S_IDP is globally identifiable (see Lemma 3 and Additional file 1 for theoretical justification).

Given a DAG G = (V,E), consider two observed nodes V _i , V _j and an unobserved node V _u . V _u may not be on any Wright’s paths between V _i and V _j . For this case,if only V _u becomes observed in O ^(k + 1), then for each observed node V _i in O ^(k), one can check whether the identifiability equation IE(V _i ,V _u ) is redundant according to Lemma 4 (see Additional file 1 for theoretical justification). That is, when none of the Wright’s paths between V _i and V _u contains detour-paths, IE(V _i ,V _u ) is redundant if and only if each Wright’s path between V _i and V _u passes at least one observed node other than V _i and V _u ; otherwise, IE(V _i ,V _u ) is redundant if and only if the Wright’s coefficient of each detour-path between V _i and V _u is globally identifiable in O ^(k) and each Wright’s path between V _i and V _u passes at least one observed node other than V _i and V _u . If V _u is on a Wright’s path between V _i and V _j , and the sufficient and necessary condition for one of the identifiability equations IE(V _i ,V _u ) and IE(V _j ,V _u ) being redundant is similar to Lemma 4 and given in Lemma 5 (see Additional file 1 for theoretical justification). Note that it can be determined whether the Wright’s coefficient of a detour-path is globally identifiable according to Lemma 2 and Lemma 3.

Based on Lemma 4 and Lemma 5, we propose a novel graphic method to calculate the identifiability gain g(V _i , O ^(k)). Let des _i denote the descendant node set of V _i , anc _i denote the ancestor node set of V _i , rel _i denote the set of nodes that are not included in des _i or anc _i . Moreover, let bound _i  ⊂ anc _i denote the boundary node set, in which every node has at least one outgoing edge to a node in rel _i . Then we can calculate g(V _i , O ^(k)) by removing the following edges from the original graph G: i) all the incoming edges to the observed nodes that are not collider nodes of detour-paths in anc _i ; ii) all the outgoing edges from some observed nodes in des _i and rel _i and these observed nodes are not the collider nodes of the detour-paths whose Wright’s coefficients are unidentifiable in O ^(k); and iii) all the outgoing edges from the observed nodes in bound _i to nodes in rel _i , and then we get a new graph denoted by G '. Let N _w denote the total number of the observed nodes that are connected with V _i via any Wright’s path in graph G '. Furthermore, one can tell from the edge-removal operation that there still exist some redundant identifiability equations in G ', because the following two types of redundancy cases have not been considered in the edge-removal operation: V _i is the downstream node of an arbitrary detour-path, and V _i is on a Wright’s path between two observed nodes in G '. Let N _r denote the number of redundant identifiability equations in G '. According to the topological structure of G ' and the node’s observation status, we can obtain N _r based on Lemma 4 and Lemma 5 (see the details in Implementation and Verification Section). It can be shown that the identifiability gain is g(V _i , O ^(k)) = N _w  − N _r (see Theorem 1 and Additional file 1 for theoretical justification).

For a given DAG G and an observation strategy O ^(k), different unobserved nodes may associate with different identifiability gains. Naturally, our strategy is to choose the unobserved node in O ^(k) with the maximum identifiability gain if it becomes observed in O ^(k + 1). However, we also recognize that, to assure that all model parameters are at least locally identifiable, certain nodes of a DAG must be observed if they are unobserved in an observation strategy (see Lemma 6 and Additional file 1 for theoretical justification). For convenience, we call such nodes the must-be-observed [14] nodes, and let \( {O}^{(0)_M} \) denote the observation strategy, in which only the MBO nodes are observed.

Lemma 1 Given a DAG G = (V,E), an observed node V _i , and an unobserved node V _u in O ^(k), if only V _u becomes observed in O ^(k + 1), the identifiability equation IE(V _i ,V _u ) is non-redundant if there exists a Wright’s path of length 1 connecting V _i and V _u .

Lemma 2 If a detour-path P has one or more exclusive upstream node, the Wright’s coefficient WP of P is globally identifiable.

Lemma 3 For a group of intersecting detour-paths, if the number of the shared upstream nodes in S_SUN is equal to or greater than the number of intersecting detour-paths in S_IDP, then the Wright’s coefficient of each detour-path in S_IDP is globally identifiable.

Then the identifiability equation IE(V _i ,V _u ) is redundant if and only if one of the above conditions holds.

Then one of the two identifiability equations IE(V _i ,V _u ) and IE(V _j ,V _u ) is redundant if and only if one of the above conditions holds.

Theorem 1 Given a DAG G = (V, E) and an unobserved node V _i in an observation strategy O, let G ' denote the sub-graph after the edge-removal operation. Then the identifiability gain is g(V _i , O) = N _w  − N _r , where N _w denotes the total number of the observed nodes that are connected with V _i via any Wright’s path in graph G ', and N _r denotes the number of redundant identifiability equations in graph G '.

Dynamic programming strategy

Let \( {O}^{(0)_G} \) denote a given observation strategy. If some of the MBO nodes are not observed in \( {O}^{(0)_G} \), \( {O}^{(0)_M} \) should be incorporated into \( {O}^{(0)_G} \) according to Lemma 6. Therefore, the initial observation strategy, denoted by O ⁽⁰⁾, should always be \( {O}^{(0)}=\left({O}^{(0)_M}\Big|{O}^{(0)_G}\right) \), where the OR operator is an element-wise operation. For example, for a DAG with six nodes, if \( {O}^{(0)_M}={\left[\begin{array}{cccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T \) and \( {O}^{(0)_G}={\left[\begin{array}{cccccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T \), then \( {O}^{(0)}={\left[\begin{array}{cccccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T \).

The dynamic programming strategy starts with the calculation of the cardinality of O ⁽⁰⁾ (that is, f(O ⁽⁰⁾)) based on Theorem 1. Specifically, let R be the number of observed nodes in O ⁽⁰⁾, V _{o_r} (r = 1, 2, ⋯, R) be the r-th observed node in O ⁽⁰⁾, and O ⁽⁰⁾{V _{o_1}, …, V _{o_r} } be the observation strategy in which only the first r observed nodes in O ⁽⁰⁾ are observed. Then \( f\left({O}^{(0)}\right)={\displaystyle \sum_{r=1}^{R-1} g\left({V}_{o\_\left( r+1\right)},{O}^{(0)}\left\{{V}_{o\_1},\dots, {V}_{o\_ r}\right\}\right)} \) can be calculated according to Theorem 1. Note that the order at which V _{o_r} is selected into O ⁽⁰⁾{V _{o_1}, …, V _{o_R} } will not change the observation strategy (e.g., O ⁽⁰⁾{V _{o_1}, V _{o_2}} = O ⁽⁰⁾{V _{o_2}, V _{o_1}}) and thus have no effect on the value of f(O ⁽⁰⁾).

The second step of our dynamic programming strategy is to define stages and their associated states. Let S denote the number of unobserved nodes in O ⁽⁰⁾, and let V _{u_s} (s = 1, 2, ⋯, S) denote the s-th unobserved node in O ⁽⁰⁾, then the dynamic programming procedure can be divided into S + 1 stages. For illustration purpose, we consider a simple example with 5 unobserved nodes, as shown in Fig. 2. The 0-th stage is actually the initialization step as described in the previous paragraph, and it has only one state, i.e., O ⁽⁰⁾. At the first stage, there are S = 5 different states; that is, only one of the unobserved nodes {V _{u_1}, V _{u_2}, ⋯, V _{u_5}} in O ⁽⁰⁾ will be selected to observe. At the second stage, since one of the five unobserved nodes has been selected at the previous stage, there are only four unobserved nodes for selection and thus four states exist (that is, {V _{u_2}, V _{u_3}, V _{u_4}, V _{u_5}}). Therefore, as shown in Fig. 2, except for stages 0 and 1, each subsequent stage has one less states than its previous stage; also, the upper triangular region (see the area above the labels of stages 1–5 in Fig. 2) is empty because the selection order of unobserved nodes does not affect the eventual observation strategy so the inclusion of such states in the upper triangular region is redundant. One can tell that the proposed stage and state definitions satisfy the optimality principle of dynamic programming [47–49].

The third step is to compute the state transition costs for searching the optimal state transition path(s). According to the definitions of stages and states, there may exist several different states at the s-th stage that can transit to the same state at the (s + 1)-th stage. For instance, four states V _{u_1}, V _{u_2}, V _{u_3} and V _{u_4} at the first stage can transit to V _{u_5} at the second stage, as shown in Fig. 2. The state transition cost from state V _{u_i} to state V _{u_j} (i ≠ j) between two consecutive stages is just the identifiability gain g(V _{u_j} , O ^(k){…, V _{u_i} , …}), where O ^(k){…, V _{u_i} , …} means that V _{u_i} is observed in O ^(k). Then the cardinality of an observation strategy can be computed by adding f(O ⁽⁰⁾) and all the state transition costs along the state transition path. Since the goal of the dynamic programming strategy is to search for the optimal observation strategies, when there exist multiple transition paths from state V _{u_i} in O ^(k) to state V _{u_j} in O ^(k + 1) (i ≠ j), the transition path associated with the maximum identifiability gain will be chosen; that is, \( f\left({O}_{V_{u\_ j}}^{\left( k+1\right)}\right)=\underset{V_{u\_ i}, i\ne j}{ \max}\left( g\left({V}_{u\_ j},{O}_{V_u{\_}_i}^{(k)}\right)+ f\left({O}_{V_{u\_ i}}^{(k)}\right)\right) \), where \( {O}_{V_{u\_ i}}^{(k)} \) is a convenient notation for O ^(k){…, V _{u_i} , …}.

It should be stressed that it is not necessary to finish all the S iterations as shown in Eq. (3). Once the cardinality \( f\left({O}_{V_{u\_ i}}^{(k)}\right) \) at the k-th stage becomes equal to or greater than the number of unknown parameters N _u , the dynamic programming process will stop and we get the SIOOR strategies.

Overview of the framework

Observation strategy design is an under-investigated problem for biological networks, despite the fact that a variety of biological networks have been actively constructed and used in numerous benchside or bedside studies. However, the existence of latent variables is a common problem due to cost, technical or other limitations, and has significantly hampered our capability to quantitatively investigate and understand such networks via, e.g., key network parameter estimation from experimental data. Identifiability analysis has long been recognized as a powerful tool to assure the accuracy and reliability of parameter estimation techniques; however, identifiability-based observation strategy design for biological networks turns out to be an unexplored field despite its substantial importance to biological network studies like structure identification.

To the best knowledge of our authors, this is the first study that tackles the problem of identifiability-based observation strategy design for biological networks described by linear SEMs. First, we introduce several new concepts such as cardinality of observation strategy and identifiability gain and mathematically formulate the identifiability-based optimal observation problem. Second, for a given network structure, the key idea is to turn a minimum number of unobserved nodes in the original observation strategy into observed such that the number of non-redundant identifiability equations becomes greater than or equal to the number of unknown model parameters (i.e., the whole system is at least locally identifiable). By counting the number of non-redundant identifiability equations, we avoid performing actual identifiability analysis on SEM and the proposed method is thus computationally efficient. Third, by defining the concepts of stage division and state transition, a dynamic programming strategy is proposed to solve the maximization problem without involving any time-consuming symbolic computation or matrix operations [33, 34]. Fourth, an efficient computing algorithm is proposed to calculate the identifiability gain of each unobserved node in a given observation strategy. More specifically, the computing process is significantly simplified by counting the number of observed nodes that connect with the node of concern via Wright’s paths after removing certain edges from the original graph.

It takes a non-constant time to compute the node identifiability gain in each iteration, and the algorithm complexity depends on the number of observed nodes. Furthermore, the number of iterations of the dynamic programming algorithm does not depend on the total number of nodes, but the number of unobserved nodes in the original observation strategy. Let S denote the number of unobserved nodes and T denote the number of observed nodes in the original observation strategy, then the computation complexity of the dynamic programing strategy is O(S ² ⋅ T).

Implementation and verification

The flowchart of the proposed algorithm for searching the structural identifiability-based optimal observation remedy is shown in Fig. 3. We have implemented the dynamic programming algorithm in R, and all the source codes and examples are freely accessible at https://github.com/Hongyu-Miao/SIOOR.

Here we describe several important technical details of the implementation. First, at the state transition step, i.e., \( f\left({O}_{V_{u\_ j}}^{\left( k+1\right)}\right)=\underset{V_{u\_ i}, i\ne j}{ \max}\left( g\left({V}_{u\_ j},{O}_{V_{u\_ i}}^{(k)}\right)+ f\left({O}_{V_{u\_ i}}^{(k)}\right)\right) \), if there exist multiple transitions that produce the same f(O ^(k + 1)), our current implementation chooses only one such transition to update the next-stage observation strategy. If needed, the R code can be slightly modified to enumerate all optimal observation strategies. Second, since the boundary node set is just a subset of the ancestor node set for a given node, the processing of the boundary nodes is incorporated into the processing of the ancestor nodes in the current implementation.

In order to verify the implementation, synthetic DAGs can be generated for this purpose, like the two DAG examples in Fig. 4. The first DAG contains 8 nodes and 13 edges, and it has only a single input node and a single output node. Moreover, the first example considers a special initial observation strategy (i.e., all nodes are unobserved) to illustrate the capability of the proposed method to design optimal observation strategy from scratch. The second DAG has multiple input and output nodes, and it considers a more general situation, that is, there exists both observed and unobserved nodes in the initial observation strategy. We analyzed the two examples using the proposed algorithm, and used the identifiability matrix method [34] to verify that the obtained observation remedies do grant (local) identifiability to all model parameters.

Applications to real biological networks

Since it is impossible to cover all the biological networks in various databases and knowledge repositories [16, 17] in one study, we choose the biological network associated with influenza A virus [39] as an application example for illustration purpose. IAV can infect birds as well as mammals including human, and it has been one of the major infectious pathogens that have caused millions of human deaths. It is thus of great scientific significance to systematically understand IAV infection and immune response mechanisms. Therefore, Matsuoka et al. [50] manually curated a comprehensive database, called FluMap, for depicting the influenza virus life cycle at the molecular level from over 500 previous publications. There are mainly five modules in FluMap: virus entry, virus replication and transcription, post-translational processing, transportation of virus proteins, and packaging and budding. Given the critical role of virus replication in influenza virus life cycle, numerous experimental studies (e.g., [42, 51–53]) have made attempts to understand virus replication mechanisms and their clinical implications. Thus, we choose to focus on the IAV replication module and analyze its observation strategy.

Since IAV replication involves many different biomolecules and complex interactions, it is usually infeasible to observe all such components and their interactions in one study. The question of concern here is how to choose a minimal number of nodes in Fig. 5 to observe such that all the model parameters become at least locally identifiable. Note that Fig. 5 is derived from Matsuoka’s work [50], and consists of 22 nodes and 26 edges; for simplicity, the catalyzers and inhibitors in this network are treated as reactants.

A relevant concept, called observability, has been previously investigated by Liu et al. [38] for complex dynamic systems. Although observability analysis also deals with observation strategies considering the existence of latent variables, it is very different from identifiability analysis in two aspects: 1) the focus of observability analysis is not model parameters but how to infer the unobserved state variables from experimentally measured outputs of a system; 2) the graphical approach proposed by Liu et al. was developed for the so-called balance equations based on mass-action kinetics, the model structures of which are very different from static linear SEMs. However, it is of interest to compare the identifiability-based observation results with those of the observability-based method. For this purpose, we assume that all the nodes in Fig. 5 are initially unobserved. After applying the proposed dynamic programing method, we get the optimal observation strategy shown in Fig. 6(a) for achieving parameter identifiability. The optimal observation strategy produced by Liu’s observability approach is shown in Fig. 6(b). According to Fig. 6(a) and (b), one can tell that the identifiability-based observation strategy contains 20 observed nodes and 2 unobserved nodes, while the observability-based strategy contains 3 observed nodes and 19 unobserved nodes. That is, for the IAV replication module, the system internal states can be inferred from a few observed output nodes if a balance equation model is used; however, it needs much more observed nodes to achieve parameter identifiability if a linear SEM is used. Such an observation is not only due to the different goals of observability and identifiability analyses, but also the differences in the underlying model structures used in observability or identifiability analyses.

Moreover, besides the nodes with an out-degree 0 or 1 as mentioned in Lemma 3, the identifiability-based observation strategy is also likely to select the nodes with a high out-degree as unobserved nodes; for instance, the two unobserved nodes viral_RNA and NP(ub) in Fig. 6(a) have the highest out-degrees 2 and 3, respectively. This is because, if an unobserved node has a high out-degree, this node is connected with many out-neighbor nodes; when its out-neighbor nodes are observed, there will exist multiple Wright’s paths that connect such out-neighbor nodes and pass this unobserved node, and the corresponding identifiability equations thus contain the parameters associated with the out-edges of this unobserved node such that these parameters can be identifiable. Interestingly, the observability-based strategy tends to select the nodes with a low out-degree as observed nodes, for example, all the nodes with 0 out-degree are observed in Fig. 6(b). It is because the nodes with an out-degree 0 in a DAG are usually the final products of chemical reactions, instead of reactants, and thus the internal states associated with other nodes can be easily inferred based on the balance equations if all the final products of chemical reactions are measured.