The partition representation of enzymatic reaction networks and its application for searching bi-stable reaction systems.

The signal transduction system, which is known as a regulatory mechanism for biochemical reaction systems in the cell, has been the subject of intensive research in recent years, and its design methods have become necessary from the viewpoint of synthetic biology. We proposed the partition representation of enzymatic reaction networks consisting of post-translational modification reactions such as phosphorylation, which is an important basic component of signal transduction systems, and attempted to find enzymatic reaction networks with bistability to demonstrate the effectiveness of the proposed representation method. The partition modifiers can be naturally introduced into the partition representation of enzymatic reaction networks when applied to search. By randomly applying the partition modifiers as appropriate, we searched for bistable and resettable enzymatic reaction networks consisting of four post-translational modification reactions. The proposed search algorithm worked well and we were able to find various bistable enzymatic reaction networks, including a typical bistable enzymatic reaction network with positive auto-feedbacks and mutually negative regulations. Since the search algorithm is divided into an evaluation function specific to the characteristics of the enzymatic reaction network to be searched and an independent algorithm part, it may be applied to search for dynamic properties such as biochemical adaptation, the ability to reset the biochemical state after responding to a stimulus, by replacing the evaluation function with one for other characteristics.

Introduction

The intracellular signal transduction system functions as a mechanism to control cell proliferation, apoptosis, differentiation, and homeostasis, and its malfunction is thought to be one of the substantial causes of cancer formation for cells. The regulatory mechanism of signal transduction systems is closely related to the mechanism of action of drugs, such as anticancer drugs, and has therefore been the subject of extensive research in recent years. In addition, from the viewpoint of synthetic biology, which has been the focus of much research in recent years, design methods for signal transduction systems have become necessary. However, because biochemical reactions are nonlinear, it is difficult to establish theoretical analysis methods, and signal transduction systems in particular have not been systematically analyzed because their interactions are more complex than those of intracellular metabolic systems, which they control. The current situation is that the parameters are fixed for each system to be analyzed, and the analysis is performed by computer simulation.

The MAPK cascade, a representative signal transduction system, which relays cell growth factor (EGF) from the cell membrane to the cell nucleus, and its abnormality is thought to be the cause of cell canceration, and much knowledge has been obtained [1-6]. A major component of the signal transduction system represented by the MAPK cascade is the cyclic reaction system, which activates and deactivates enzymes by phosphorylation and dephosphorylation. The cyclic reaction system is a combination of two post-translational modification reactions. Therefore, we have been conducting a comprehensive analysis of the stability of regulatory networks that consist of activation and inactivation cyclic reaction systems of enzymes and their mutual controls [7]. There, the control relationship between the cyclic reaction systems is represented by a control matrix. The control matrix, which is an adjacency matrix, has 1 as its i-th row and j-th column components when the activating enzyme of cyclic reaction system j catalyzes the activation of cycle reaction system i, and -1 when it catalyzes the deactivation.

In a similar study, Kuwahara et al. [8] exhaustively analyzed the influence of the control structure on the stochastic properties for regulatory networks of three to five nodes with only one feedback control. They also mentioned the control structures in which bistability emerges. Ramakrishnan et al. [9], Shah et al. [10], and Siegal-Gaskins et al. [11] have conducted exhaustive analyses of the ultrasensitive properties and bistability of signaling systems. Ultrasensitivity was reported as a response characteristic of the MAPK cascade, a typical signaling system mentioned above [1]. In a system with ultrasensitivity, the steady-state enzyme concentration of the output chemical species with respect to the enzyme concentration of the input chemical species increases rapidly after a certain threshold of the input. Bistability is when this ultrasensitivity is more strongly expressed, and hysteresis appears in the change of output to input. Thus, in a region with input, two steady state values of the output appear. In particular, it is said to be resettable if it has the property of being able to mutually transition between two steady states by varying the total concentration of the input chemical species within the range of protein concentrations expected in the cell [12].

Ma et al. [13] and Yao et al. [14] analyzed the biochemical adaptation and robustness properties using a similar control network of three nodes. Biochemical adaptation is the property that when the concentration of an input chemical species is increased in a steady-state system, the response of the output chemical species transiently increases and then decreases to a value close to the original steady-state value [13]. Robustness is the property that the various response properties of a system are less sensitive to the concentration of other chemical species in the system [12].

Adler et al. [15] focused on fold-change detection, which is a dynamic property in which the temporal pattern of the system’s output in response to a transient input change in the cell depends on the ratio rather than the difference of the input changes, demonstrating the effectiveness of comprehensive analysis using regulatory networks. In these studies, the enzymatic reaction mechanism is mainly based on the Michaelis-Menten approximation or a simplified linear equation in order to reduce the computational complexity. However, it has been reported that approximations can lead to inaccuracies in the bistability and dynamics of the system [16, 17]. In order to accurately analyze the properties of a system, it would be better to describe the system using only the law of mass action, but this would increase the number of parameters of the system and make the originally huge search space even larger, making exhaustive analysis difficult to carry out. Therefore, there is a need to develop a method to search for enzymatic reaction networks with useful properties such as ultrasensitivity, bistability, biochemical adaptation, robustness, and fold-change detection as mentioned above [18-20].

In this study, we focus on the enzymatic reaction network, which is formulated as a regulatory network in which post-translational modification reactions are the elements and they mutually regulate each other. Post-translational modification reactions, such as phosphorylation, are smaller units than the cyclic reaction systems. Instead of the control matrix described above, we propose to use a partition of the set as its representation, which is considered to be more suitable for search. In addition, we report the results of our search for a bistable and resettable enzymatic reaction network with the aim of demonstrating its effectiveness. Resettable is the property of being able to transition between two steady states by changing the total concentration of the input chemical species within the range of expected protein concentrations in the cell [12].

Materials and methods

Partition representations for enzymatic reaction networks

A post-translational modification reaction with Michaelis-Menten-type enzymatic reaction as the reaction mechanism is considered as the basic component, considering to construct the enzymatic reaction network composed of N components. In this regard, reactions in which multiple enzymes bind to a single substrate simultaneously are not considered for simplicity in this formulation. The Michaelis-Menten-type enzymatic reaction is a series of chemical reactions in which an enzyme (E) binds to a substrate (S) to form a temporary enzyme-substrate complex (C), the substrate is transformed into a product (P) on the complex, and then decomposes into the product (P) and the enzyme (E), as shown in Eq (1), where the reaction rate constants are denoted by a, d, and k. The substrates, enzymes, and products of N Michaelis-Menten-type enzymatic reactions are denoted by Si, Ei, and Pi, respectively, and the set is denoted by M. That is, M = {S1, E1, P1, S2, E2, P2, …, SN, EN, PN}.

Take one partition P of the set M. A partition is a family of mutually disjoint subsets of a set, and the union of those subsets is the original set. By identifying the elements of the partition, i.e., the chemical species that belong to one subset of the original set M, an enzymatic reaction network is constructed from multiple post-translational modification reactions. Fig 1 shows an example of a basic enzymatic reaction network consisting of post-translational modification reactions. The primary building blocks of post-translational modification reactions are numbered consecutively; the substrate, enzyme, and product of the i-th post-translational modification reaction are labeled Si, Ei, and Pi, respectively. If the chemical species to be equated is the substrate Si and the product Pi, the names of the species to be equated are listed in the rectangle. The enzyme Ei is represented by the small circle between the red arrows, and if there is a chemical species that is identical to the enzyme, it is connected to that chemical species by a dotted line. It may be easier to understand intuitively if you think of the small circle as representing the enzyme-substrate complex.

Fig 1

Basic enzymatic reaction networks consisting of post-translational modification reactions.

Examples of enzymatic reaction networks composed of post-translational modification reactions as the primary building blocks; a: product of enzymatic reaction 1 catalyzes enzymatic reaction 2, b: two-step reaction, c: branching reaction, d: confluence reaction, e: autocatalysis, f: association-dissociation reaction, g: dimer formation separation reaction. The primary building blocks of post-translational modification reactions are numbered consecutively; the substrate, enzyme, and product of the i-th post-translational modification reaction are labeled Si, Ei, and Pi, respectively. If the chemical species to be equated is the substrate Si and the product Pi, the names of the species to be equated are listed in the rectangle. The enzyme Ei is represented by the small circle between the red arrows, and if there is a chemical species that identifies with the enzyme, it is connected to that species by a dotted line.

Fig 1A–1D are examples of an enzymatic reaction network consisting of N = 2, i.e., two post-translational modification reactions, where the original set of partitions is M = {S1, E1, P1, S2, E2, P2}. Fig 1A shows the enzymatic reaction network represented by the partition {{S1}, {E1}, {P1, E2}, {S2}, {P2}} of M, where the product of enzymatic reaction 1 functions as the enzyme in enzymatic reaction 2. The identification of Pi and Ej corresponds to the fact that the product of enzymatic reaction i functions as the enzyme of enzymatic reaction j. This is a typical mode of enzymatic reaction chain in which the activated enzyme catalyzes other enzymatic reactions as seen in signal transduction systems. Fig 1B shows the enzymatic reaction network represented by the partition {{S1}, {E1}, {P1, S2}, {E2}, {P2}} of M, where the products of enzymatic reaction 1 are the substrates of enzymatic reaction 2. The identification of Pi and Sj is a representation of the two-step enzymatic reaction from Si to Pj. Fig 1C and 1D show the enzymatic reaction networks represented by the partitions, {{S1, S2}, {E1}, {P1}, {E2}, {P2}} and {{S1}, {E1}, {S2}, {E2}, {P1, P2}} of M, respectively. The identification of Si and Sj, and Pi and Pj, respectively, is a representation of the branching into and confluence from two enzymatic reactions.

Fig 1E to g are examples of enzymatic reaction networks consisting of N = 1, i.e., one post-translational modification reaction, where the original set of partitions is M = {S1, E1, P1}. Fig 1E shows the enzymatic reaction network represented by this partition {{S1}, {E1, P1}} of M, which is a representation of an autocatalytic reaction. It corresponds to the case where i = j = 1 in Fig 1A. Fig 1F shows the enzymatic reaction network represented by the partition {{S1, P1}, {E1}} of M, showing the association and dissociation reaction between the identical substrate and product S1 = P1 and the enzyme E1. The identification of Si and Pi is a representation of the association and dissociation reaction with the enzyme. Fig 1G shows the enzymatic reaction network represented by the partition {{S1, P1, E1}} of M, where the identified substrate and product and the enzyme S1 = P1 = E1 associate and dissociate themselves, corresponding to the dimer formation separation reaction. The identification of Si, Pi, and Ei is an expression of the dimer formation-separation reaction.

Fig 2 shows examples of more complex enzymatic reaction networks and their partition representations. Fig 2A shows an example of an enzymatic reaction network consisting of N = 4, i.e., four post-translational modification reactions, where the original set of the partition is M = {S1, E1, P1, S2, E2, P2, S3, E3, P3, S4, E4, P4} and the partitions is {{S1, P2}, {P1, S2, E1, E4}, {S3, P4}, {P3, S4, E2, E3}}. This is an example of an enzymatic reaction network consisting of two cyclic reaction systems with auto-activating and mutually inhibitory feedbacks, which is an typical of enzymatic reaction network with bistable properties. The MAPK cascade, one of the representative signaling systems, is shown in Fig 2B. The MAPK cascade consists of a cascade of three kinases, MAPKKK, MAPKKK, and MAPK. The activated enzyme activates the phosphotransferase in the next step, thereby transmitting signals within the cell. An example of an enzymatic reaction network consisting of 10 post-translational modification reactions, with M = {S1, E1, P1, S2, E2, P2, …, S10, E10, P10} can be represented by the partition {{P2, S1}, {P1, S2, E3, E4}, {E1}, {E2}, {P6, S3}, {P3, P5, S4, S6}, {P4, S5, E7, E8}, {P10, S7}, {P7, P9, S8, S10}, {P8, S9}, {E5}, {E6}, {E9}, {E10}}, where {P2, S1} and {P1, S2, E3, E4} are the inactive and active MAPKKKs, respectively. Also, {P6, S3}, {P3, P5, S4, S6}, and {P4, S5, E7, E8} are inactive, one phosphorylated, and two phosphorylated MAPKKs, respectively. {P10, S7}, {P7, P9, S8, S10}, and {P8, S9} are inactive, one-phosphorylated, and two-phosphorylated MAPKs, respectively.

Fig 2

Complex enzymatic reaction networks consisting of post-translational modification reactions.

The way to represent the figure is the same as in Fig 1; a: auto-activating mutual inhibitiory network, b: MAPK cascade.

The advantage of the partition representation is that every partition corresponds to an enzymatic reaction network, and thus there is a one-to-one relationship between the partition and the enzymatic reaction network. The number of different enzymatic reaction networks consisting of N post-translational modification reactions is the same as the total number of different partitions of the set with number of elements 3N, a number whose value is known as the Bell number B, expressed in the recurrence equation as in Eq (2), where B0 = B1 = 1. For example, the total number of enzymatic reaction networks consisting of four post-translational modification reactions, as shown in Fig 2A, is B12 = 4213597.

In the space where the partition is an element, we can introduce the distance D (A, B) between the partition A = {a1, a2, …, am} and the partition B = {b1, b2, …, bn} of the set M as shown in Eq (3). Here, ai and bj are subsets of the set M. The symbol # is a function that returns the number of elements in a set. Therefore, # (a ∪ b—a ∩ b) is the number of elements that are not common to the sets a and b. For example, let the source set of the partition be M = {x, x, x, x, x, x}, and consider three partitions, A = {a, a, a}, B = {b, b}, and C = {c, c}, where a = {x, x, x}, a = {x, x}, a = {x}, b = {x, x}, b = {x, x, x, x}, c = {x, x, x}, c = {x, x, x}. Then d (A, b) = min s (a, b) = s (a, b) = 1, d (A, b) = min s (a, b) = s (a, b) = 2, so max d (A, b) = d (A, b) = 2, and similarly max d (B, a) = 3. Therefore, D (A, B) = max (2, 3) = 3. Similarly, we can see that D (B, C) = 1 and D (A, C) = 2. In other words, A and B are the farthest, B and C are the closest, and A and C are in the middle distance. To check this in an enzymatic reaction network, for example, let the original set of chemical species be M = {S1, P2, E2, P1, E1, S2}. In other words, x = S1, x = P2, x = E2, x = P1, x = E1, x = S2. Then partition A becomes {{S1, P2, E2}, {P1, E1}, {S2}}, B becomes {{S1, P2}, {E2, P1, E1, S2}}, and C becomes {{S1, P2, E2}, {P1, E1, S2}}, corresponding to a, b, and c of the enzymatic reaction network shown in Fig 3.

Fig 3

Enzymatic reaction networks specified by divisions and distances between the divisions.

The way to represent the figure is the same as in Fig 1. The enzymatic reaction networks represented by the partitions; A: {{S1, P2, E2}, {P1, E1}, {S2}}, B: {{S1, P2}, {E2, P1, E1, S2}}, C: {{S1, P2, E2}, {P1, E1, S2}}.

It can be seen that the similarity of the enzyme reaction network correlates with this distance. Furthermore, this distance between partitions satisfies the triangle inequality and can be used to visualize the distribution of partitions.

Derivation of differential equation systems and conservation laws from partitions

From the partition, we can derive a system of differential equations and conservation laws that describe the behavior of the corresponding enzymatic reaction network as follows. If the number of post-translational modification reactions constituting the enzymatic reaction network is N, the total number of elements in the set from which the partitioning is made is 3N, since three elements, that is, a substrate, an enzyme, and a product, correspond to one post-translational modification reaction. Therefore, if we denote the elements by x, the source set of the partition can be expressed as M = {x, x, …, x}.

Let one of the partitions of M be A = {a, a, …, a}, and let s be the set of subscript i of x which is an element of the kth element a of the partitions A. For this partition A, a system of differential equations can be generated by the following procedure.

Generate 4N reaction rate equations for substrate Si, enzyme Ei, enzyme-substrate complex Ci, and product Pi in each post-translational modification reaction, as shown in Eq (4). Note that x, x, …, x correspond to the enzyme-substrate complex Ci.

For each element ak in the partition, generate a chemical species yk that identifies the all elements of ak, and replace the original reaction rate equations for the elements of a with the following Eq (5), while all elements of ak appearing on the right side of Eq (5) are replaced by yk. The s is the subscript set of the elements of a.

For the derivation of Eq (4), if only the law of mass action is used as the reaction equation for the post-translational modification reaction, and the variable name is expressed as x[chemical species], the reaction equation derived is the following for the post-translational modification reaction i. For example, if the original set of partitions is M = {S1, E1, P1, S2, E2, P2}, and the partition is A = {{S1}, {E1}, {P1, E2}, {S2}, {P2}} as seen in Fig 1A, then by applying the law of mass action, the above Eq (4) becomes the following Eq (7), where the variable name is represented by x[chemical species]. The equations corresponding to Eq (5), obtained by equating P1 and E2 according to division A, are Eq (8). In particular, when some chemical species are considered identical, they are listed in brackets. For the conservation law, we first generate a set of lists Q = {q1, …, qi, …, q2} of chemical species whose total concentration is conserved for each post-translational modification reaction. The reason why the total number is 2N is that for each post-translational modification reaction, there is a conservation law for the enzyme and a conservation law for the substrate. Next, the following procedure is applied sequentially to each element ak of the partition A, in turn.

Concatenate all the list of elements in Q that contain elements in common with a, and let r be the list.

Remove the elements of a from the list rk and add the chemical species y, which is identical to the elements of a.

Add r to the element of Q that does not contain an element in common with a, and make it Q again.

It is important to note that the above procedure concatenates, rather than sums, the elements of Q that contain elements in common with ak. This is because the overlap is intrinsic when the list to be combined contains the common enzyme-substrate complex Ci.

If only the law of mass action is used as the reaction equation for the post-translational modification reaction, and the variable name is represented by x[chemical species], the initial elements of the set of conservation laws Q to be derived are the following for the post-translational modification reaction i. For example, for the enzymatic reaction networks in Fig 1A, the partition is A = {{S1}, {E1}, {P1, E2}, {S2}, {P2}}. Then, if the variable names are denoted by x[chemical species] as in the example above, the original list of conservation laws becomes {{x[S1], x[C1], x[P1]}, {x[E1], x[C1]}, {x[S2], x[C2], x[P2]}, {x[E2], x[C2]}}, and by equating P1 and E2, we get Q = {{x[S1], x[C1], x[P1, E2], x[C2]}, {x[E1], x[C1]}, {x[S2], x[C2], x[P2]}}.

Resettabe bistability

A search for bistable and resettable enzymatic reaction networks was attempted to demonstrate the effectiveness of the partition representation of enzymatic reaction networks. Resettable is the property of being able to mutually transition between two steady states by varying the total concentration of the input chemical species within the range of expected protein concentrations in the cell [12]. A typical resettable bistable enzymatic reaction network is shown in Fig 2A.

The relationship between bistability and resettability when the input chemical species is E1 and the output chemical species is P2 is shown in Fig 4. For example, in the enzymatic reaction network of Fig 2A, the input chemical species E1 corresponds to the small circle labeled 1, and the output chemical species P2 corresponds to the rectangle in the lower left corner that equates S1 and P2. The horizontal axis is the logarithm of the input species concentration with a base of 2, and the vertical axis is the total output species concentration normalized to between 0 and 1. The range of protein concentrations expected in the cell is 2−5 to 25 in m mol/m3.

Fig 4

Relationship between bistability and resettability.

The horizontal axis is the total concentration of the enzyme including the input chemical species E1 and the vertical axis is the relative concentration of the output chemical species P2 to the total concentration; a: resettable bistable, b: non-resettable bistable, but the transition from the upper stable state to the lower can occur, c: non-resettable bistable. The solid line corresponds to the stable equilibrium point, i.e., the steady state value, and the dashed line corresponds to the unstable equilibrium point, The solid blue lines at both ends depict the monostable state.

The curve represents the concentration of the output species at equilibrium, where the solid line corresponds to the stable equilibrium point, or steady state value, and the dashed line corresponds to the unstable equilibrium point.

Fig 4A shows the resettable bistable aspect. The solid blue curves at both ends of the graph correspond to a single steady state, which are monostable states. The red curve in the center has two steady states, corresponding to bistable states. As the total concentration of the input species is increased from the smaller one, the black arrow on the right makes a discontinuous jump from the upper steady state to the blue monostable state on the right. When the total input species concentration is reduced from that state, the system jumps to the steady state below at the black arrow on the left side. Thus, it is a resettable bistable because it can mutually transition between two steady states.

Fig 4B is an example where the curve in Fig 4A extends to the right and the region of monostability on the right exceeds the range of protein concentrations expected in the cell. In this case, the black arrows show that it is possible to jump from the upper steady state to the lower steady state, but it is not possible to jump from the lower steady state to the upper steady state, even if the total input chemical species concentration is increased within the range of expected protein concentrations in the cell. This is an example of a bistable that is not resettable.

Fig 4C shows an example where the graph extends further to the left. In such a case, if the initial state of the system is above the red dotted line in the center, it shifts to the upper steady state, and if it is below, it shifts to the lower steady state, and the state cannot be changed even if the total concentration of the input chemical species is changed within the range of the expected protein concentration in the cell. This is another example of bistability that is not resettable.

Exploration of enzymatic reaction networks in the partition representation space

For the partition representation of the enzymatic reaction network, an algorithm to search for the enzymatic reaction network with the desired response characteristics can be implemented in a natural way by sequentially applying the following two partition modifiers.

Separate partition modifier: Randomly selects an element from a partition and randomly separates it into two elements. For example, by separating the second element {P1, E1, S2} of the partition {{S1, P2, E2}, {P1, E1, S2}} corresponding to the enzymatic reaction network C shown in Fig 3 into two elements {P1, E1} and {S2}, we obtain the partition {{S1, P2, E2}, {P1, E1}, {S2}} corresponding to the enzymatic reaction network A is obtained.

Join partition modifier: Combines two randomly selected elements. For example, in the example of Fig 3 above, by combining the second and third elements of the partition corresponding to the enzymatic reaction network A, the partition corresponding to the enzymatic reaction network C is obtained.

First, we set up an evaluation function ϕ(P,ν) that quantifies the desired response properties of the dynamics or steady state of the enzymatic reaction network generated from the partition. ϕ(P,ν) is a function of the partition P, with the parameter ν, which is a pair of reaction rate constants of each post-translational modification reaction and the total concentration of enzymes that constitute the enzymatic reaction network corresponding to the partition.

In the following algorithm, we use the function Ф(P) to perform a random search on the parameter value ν for a given partition P. The function Ф(P) quantifies the evaluation value at λ randomly chosen parameter values from the set of candidate values Λ of the parameter ν by ϕ(P,ν) and returns the pair of the largest evaluation value and the parameter ν at that time as the value.

The main loop of the algorithm is as follows.

Determine one partition P at random.

If the value of the evaluation function Ф(P) is less than ρ or the number of iterations is less than γ, repeat P = ψ(P, P, σ).

Return P as the value.

The partition search function ψ(P, P, σ) is a recursive function whose arguments are the initial partition P, the partition P, and the partition depth σ described below. ψ(P, P, σ) returns the partition whose evaluation function is greater than or equal to Ф(P).

If σ is zero, return P as the value.

Generate P’ from P with the join partition modifier.

If Ф(P’) > Ф(P), return P’ as the value.

Generate P” from P with the separate partition modifier.

If Ф(P”) > Ф(P), return P” as the value.

let P = ψ(P, P’, σ-1).

If Ф(P) > Ф(P), return P as the value.

let P = ψ(P, P”, σ-1).

If Ф(P) > Ф(P), return P as the value.

If none of the above conditions are satisfied, return P as the value.

We apply this random search algorithm to the aforementioned search for enzymatic reaction networks with resettable bistability with the aim of demonstrating the effectiveness of the partition representation of enzymaric reaction networks. In other words, we design an evaluation function ϕ(P,ν) for resettable bistability and use it in this algorithm.

Results

Evaluation function for resettable bistability

Bistability is the property of having two steady state values for a given parameter value ν, as described above. Furthermore, resettable is the property that allows the transition between these two steady states by changing the value of the total concentration of the enzyme containing the input chemical species. Here, we create an evaluation function ϕ(P,ν) that gives a high evaluation value for a typical resettable bistable aspect as shown in Fig 4A, i.e., an enzymatic reaction network that is monostable at the edge of the range taken by the total concentration of input chemical species and bistable in the central part.

ϕ(P,ν) is a function of the partition P, with the reaction rate constants of each post-translational modification reaction and the total concentration of enzymes comprising the enzymatic reaction network corresponding to the partition as parameters ν. The variations in the total concentration of enzymes containing the input chemical species were set to 21 discrete values of Ω = 2−5, 2−4.5, …, 25. These values, as well as the reaction rate constants to be set as ν, are set to include approximately the values reported for the kinases that make up the MAPK cascade, which is known to be a typical cellular signaling system.

The following procedure corresponds to the aforementioned evaluation function ϕ(P,ν). P and ν are given as arguments.

For a given parameter value ν, only the total concentration of the enzyme containing the input chemical species is varied sequentially and discretely for 21 values of Ω = 2−5, 2−4.5, …, 25, the total concentration of each enzyme is randomly assigned as the initial value so as to satisfy the conservation law of each enzyme. The steady state values were obtained by numerically solving the generated differential equations and making convergence judgments as appropriate. The details are described in the next section. Then 21 pairs of steady state values of the output chemical species are obtained. Here, the steady-state value of the output enzyme species is normalized to its total concentration, and is a value between 0 and 1.

This is repeated θ times to obtain θ pairs of normalized steady state values for each of the 21 different total concentrations of the input chemical species, where only data for which all 21 values are converged are used. Therefore, the number of them is at most θ. Fixing one from 21 different total concentrations of the input chemical species means that we have at most θ steady state values. When all these values are the same, it means monostable, and when they are divided into two different values, it means bistable.

The degree of bistability can be quantitatively evaluated by the variance of several steady-state values obtained from different initial states. This value takes a maximum value of 1/4 in a perfect bistable state, i.e., when half of the values are 0 and the rest are 1, and a minimum value of 0 in a monostable state, i.e., when all the values are the same. It also takes an intermediate value in the mixed states. By multiplying the value of the variance by four, we can normalize the number representing the degree of bistability from 0 to 1.

Up to this point, we obtain 21 normalized variances V = {V, V,…,V}, which indicate the degree of bistability for each of 21 values of the total concentration of the enzyme containing the input chemical species: 2−5, 2−4.5, …, 25.

Return the value of the function eVL for this V, eVL(V), as the evaluation value.

The function eVL is the function defined in Eq (10). Resettable bistability is the property of being monostable at the edge of the range taken by the total concentration of the input chemical species and bistable in the middle part, as explained in Fig 4. The eVL evaluates the degree of this property. Its input V is a vector of the above 21 quadruple variance values and Vi denotes the each element. The value of the function eVL was normalized by a maximum value so that it would be between 0 and 1. The α is the coefficient for this purpose, which is the reciprocal of the maximum value 0.604762 that the summation equation on the right side takes when V = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0}. The 3D shape of the function wVL that makes up the function eVL is shown in Fig 5. Its input v is the quadruple variance value V and l is the position i of V in the vector V. The function wVL is designed to be maximal for perfect monostability (v = 0) at the edge of the region, and also to be maximal for perfect bistability (v = 1) at the center of the region, and continuous intermediate evaluation values in the middle of the region.

Fig 5

Function wVL for evaluation of resettable bistability.

Input v is the quadruple variance value V and l is the position i of V in the vector V. The wVL is designed to take a maximum value for perfect monostability (v = 0) at the edge of the region, a maximum value for perfect bistability (v = 1) in the center, and continuous intermediate evaluation values in the middle of the region.

Hyperparameters of the search algorithm

The following is a summary of the hyperparameters included in the algorithm for searching for resettable bistable enzymatic reaction networks, which is an application of the partition representation, and the values used in the actual search. First, the hyperparameters that are independent of the characteristics of the enzymatic reaction networks to be explored include the following.

Number of post-translational modification reactions N: 4

Input and output chemical species: input is E1 and output is P2

Partitioning depth (depth of recursive search in binary tree) σ: 4

Maximum number of loops to change the division γ: 150

Value of the evaluation function to be judged as convergence ρ: 0.8

A set of candidate values of a parameter to try randomly Λ: {2−5, 2−4.5, …, 25}

Number of pairs of parameter values to try randomly λ: 5

As mentioned earlier, the enzymatic reaction network shown in Fig 2B is known as a typical resettable bistable enzymatic reaction network [7], and the number of post-translational modification reactions that make up the network is 4, so the value of N was set to 4. The input and output chemical species are set to be E1 and P2, respectively, which belongs to the different post-translational modification reactions, because the enzymatic reaction network explored would be severely limited if the input and output chemical species were those of the same post-translational modification reaction.

We did not know whether the values we chose for the depth of recursive search σ and the maximum number of loops γ were optimal or not. After trying several combinations, we chose the one that would allow us to explore a single enzyme reaction network in at least a few days. The set of values Λ used for the total concentration of enzymes and the rate constants of the enzymatic reactions was set to include approximately the values in mmol/m3-s system reported for the phosphatases comprising the MAPK cascade [21-25]. For λ, we tried several other values, but did not find much difference.

Hyperparameters specific to the search target, i.e., the search for bistable enzyme reaction networks, include the following. In general, these will alter for different search targets.

A series of total concentrations of the input chemical species Ω: {2−5, 2−4.5, …, 25}

Number of steady-state values with random initial values θ: 5

We used the same variation of Λ as described above for Ω. We tried several other values, but did not find much difference for θ.

To obtain the above steady state values, the generated differential equations were solved numerically many times while updating the initial values, and the convergence was judged. This part has the following hyperparameters, which are common when steady state values are needed in the search process.

Calculation time for one step τ: 10 s

Maximum number of step repetitions κ: 360

Convergence decision error ε: 10−3

Wolfram Mathematica’s NDSolve function was used to solve the differential equations numerically. The protein concentration after τ seconds was calculated by NDSolve, and was judged to be converged when all of the relative error of the change, i.e., the ratio of the change to the total concentration of the protein, was less than ε. If there is no convergence, the protein concentration after another τ seconds is calculated and the same decision is made. If it does not converge after repeating this process κ times, it returns the result that it did not converge. The product of τ and κ, 3600 seconds, or 1 hour, was taken as the critical time to reach steady state. This time was taken as a reference for the response time of the signaling system in the cell. The values of the calculation time τ for one step and the convergence decision error ε were adopted after trial and error.

Enzymatic reaction networks found by the search algorithm

The search was conducted on a machine with a 3.4GHz i7-3770 processor CPU, 16GB of memory, and a Windows 10 operating system. The search program was implemented in Wolfram Mathematica 12.0. The bistable enzymatic reaction networks found in the search at the values of the hyperparameters described above are shown in Fig 6A–6H. Each of them uses a different seed of random numbers for the initial setting. The maximum CPU time required to discover a single enzymatic reaction network was approximately 30 hours. For example, the enzymatic reaction network in Fig 6A consists of three chemical species: the input is the chemical species that identifies E1, E2, P1, P4, S1, and S4, and the output is the chemical species that identifies P2 and S3. The rest is a chemical species that identify E3, E4, P3, and S2. In particular, the input chemical species is in equilibrium with its dimer. Fig 6G is identical to the typical bistable enzymatic reaction network shown in Fig 2A. The reaction rate constants for each enzymatic reaction network, the total concentration of each chemical species, and the derived differential equation systems are shown in S1–S4 Tables provided as the supporting information, respectively.

Fig 6

Bistable enzymatic reaction networks found in the search.

The enzymatic reaction networks a—h use the different seeds of random numbers. The way to represent the figure is the same as in Fig 1.

Fig 7 shows the bistable aspect of each enzymatic reaction network, where a—h corresponds to a—h in Fig 6. In order to plot the bistable aspect, the Monte Carlo method was used to find the steady state values of the output chemical species for each value of the total concentration of the enzyme containing the input chemical species. That is, when numerically solving the system of differential equations derived from the partition, the initial values were randomly allocated to the chemical species included in each conservation law. We can see that all the obtained enzymatic reaction networks exhibit resettable bistability. Resettable is the property of being able to go back and forth between two bistable steady state values by changing the value of the input.

Fig 7

Bistable aspects of enzymatic reaction networks found in the search.

The enzymatic reactoin networks a—h correspond to Fig 6. The horizontal axis is the total concentration of the enzyme including the input chemical species, E1, in mmol/m3, and the vertical axis is the relative concentration of the output chemical species, P2, to the total concentration.

Fig 8 shows an aspect of the search process. Fig 8A shows the evolution of the evaluation value as the search progresses. The horizontal axis is the number of times the structure or parameter values of the enzymatic reaction network were changed. The vertical axis is the evaluation value, with the maximum being 1. The legends a—h corresponds to a—h in Fig 6. Since the evaluation value reaches 0.6 in all cases in the first few steps, and almost all graphs overlap up to that point, the range of the vertical axis is set to 0.55 or higher so that there is less overlap thereafter. The search is limited to 150 times, and the search is terminated when the evaluation value exceeds 0.8. It can be seen that the aspect of convergence varies depending on the seeds of random numbers used at initial settings. Fig 6B shows movements within the partition space, based on the distance between the partitions defined by Eq (3). Wolfram Mathematica’s Graph function was used for drawing. The distance between all partitions was measured by the introduced distance function and the values were used in the EdgeWeight option. The GrpahLayout option specifies SpringEmbedding, which is arranged so that the total energy of the mutually coupled partitions is minimized by a spring of force equivalent to the distance between the partitions. The legends a—h corresponds to a—h in Fig 6. The red circles denote the final partitions in each search. We can see that the search is moving through the entire partition space.

Fig 8

Aspects of the search process.

The enzymatic reaction networks a—h correspond to Fig 6. A: Aspect of convergence. Horizontal axis is the number of iterations of the search. The vertical axis is the evaluation value. Convergence occurs when the evaluation value is greater than 0.8 or the number of iterations exceeds 150. B: Distribution of partitions in the search process based on the distance between partitions. The search proceeds in the direction of the arrow, and the red circles are the final points.

Discussion

We have shown that an enzymatic reaction network composed of post-translational modification reactions can be represented by a partition of the set whose elements are the enzymes, substrates, and products. By identifying elements of each subset in the partition as the same chemical species, an enzymatic reaction network can be constructed. In particular, by equating the substrate and product within the same post-translational modification reaction, we can also describe the association and dissociation reaction between enzyme and substrate (product). However, it does not represent the binding of multiple proteins as found in scaffold proteins and as seen in membrane receptor proteins. As an extension of the partition representation, enzyme-substrate complexes may be added as elements of the original set, but an algorithm for deriving the system of differential equations and conservation laws from the partition needs to be devised.

Furthermore, the partition representation of enzymatic reaction networks was applied to the search for bistable networks. The two partition modifiers introduced into the partition representation worked well and various bistable enzymatic reaction networks were obtained. One of the enzymatic reaction networks obtained in the search were the typical bistable enzymatic reaction network in which two cyclic reaction systems with positive auto-regulatory feedbacks mutually negatively regulate each other.

Regarding the performance of the search algorithm, it was confirmed that the search was completed within an acceptable processing time. It is especially noteworthy that the processing time required for the search was greatly reduced compared to an exhaustive search. In addition, the visualization of the search path using the distance between the partitions introduced in this study was able to show the aspect of the search.

It is not difficult to increase the granularity of the primary building blocks from post-translational modification reactions to cyclic reaction systems. Instead of the original set of elements being the substrate, enzyme, and product of the post-translational modification reaction, we can have four types of enzymes: the active enzyme, the inactive enzyme, the activating enzyme that catalyzes the reaction that turns the inactive enzyme into the active enzyme, and vice versa. We plan to try it in the future.

For the hyperparameter of the search algorithm, which is the search depth, and the parameter values for evaluating the structure of an enzymatic reaction network, which are the reaction rate constant and the number of times the total concentration of each chemical species is randomly tested, we set the search depth to 4 and the number of random tests to 5, but we did not systematically test other values. It is possible that adjusting these values can speed up the process. Further speed-up can be achieved by using GPUs or by using tQCCM [16], an improved form of the Michaelis-Menten approximation, as a formulation of the reaction mechanism instead of the law of mass action.

Although the bistability targeted in this study is a steady-state property, the search algorithm is divided into an evaluation function specific to the characteristics of the enzymatic reaction network to be explored and an independent algorithmic part, so that it may be possible to apply it to the exploration of dynamic properties such as biochemical adaptation, by replacing the evaluation function. We plan to try it in the future. The number of post-translational modification reactions, which is the primary building block, was set to be 4, but even when the number was set to be 3, a bistable enzymatic reaction networks could be found. If the mechanism for automatic increase or decrease in the number of post-translational modification reactions can be included in the search algorithm, it will be possible to find a more optimal enzyme reaction network.

An example of the transitions in search for enzymatic reaction networks.

Transitions of the enzymatic reaction networks in the process of searching for the discovered enzymatic reaction network Fig 6G. The top left is the initial enzymatic reaction network and transitions to the right, then moves to the second line and transitions from left to right, and the third line transitions from left to right as well. The bottom right is the final enzymatic reaction network Fig 6G.

(EPS)

Click here for additional data file.

Reaction rate constant values for the enzymatic reaction networks found.

Each column is the value of the respective reaction rate constant. Each line corresponds to a—h in Fig 6.

(XLSX)

Click here for additional data file.

The total concentration of enzymes in each enzymatic reaction network found.

Each line corresponds to a—h in Fig 6. The second column is a list of sets of chemical species corresponding to the conservation laws. The third column is their total concentration in the same order as the second column. However, for the set containing E1, the input chemical species, we have moved it in the range of 2-5 to 25, instead of the values in the table.

(XLSX)

Click here for additional data file.

Differential equation system for each enzymatic reaction network found (a—d).

The first column corresponds to a—d in Fig 6. The second column is the name of the variable to be differentiated on the left side of the differential equation and the third column is the right side of the differential equation. The variable name is represented by x[chemical species name]. In particular, when some chemical species are considered identical, they are listed in brackets.

(XLSX)

Click here for additional data file.

Differential equation system for each enzymatic reaction network found (e–h).

The first column corresponds to e–h in Fig 6. The second column is the name of the variable to be differentiated on the left side of the differential equation and the third column is the right side of the differential equation. The variable name is represented by x[chemical species name]. In particular, when some chemical species are considered identical, they are listed in brackets.

(XLSX)

Click here for additional data file.

(ZIP)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

15 May 2021

PONE-D-21-10250

The partition representation of enzymatic reaction networks and its application for searching bistable reaction systems

PLOS ONE

Dear Dr. Naka,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Please note that both referee find the manuscript hard to follow and identify unsupported statements. We expect that the presentation of this work should be significantly improved before we can contact the referees again.

Please submit your revised manuscript by Jun 29 2021 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see:  http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols . Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at  https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols .

We look forward to receiving your revised manuscript.

Kind regards,

Ivan Kryven

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: No

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: No

Reviewer #2: N/A

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The paper is not intelligible the way it is currently written. It is difficult to say whether the findings merit publication if written better, because in the current form it is difficult to understand even what the main point is. The author says "proposed search algorithm worked well" in the abstract without specifying the basis for this statement. Similarly, there is a vague and general remark at the end of the abstract "seems to be applicable to the search for dynamic properties".

But mostly the main text of the paper is extremely unclear. Maybe the author can summarize the goal of the paper in a short paragraph, and explain one example clearly. Unfortunately, the way it is currently written, I am afraid the paper may be of little use to anyone who is not completely well-versed with the author's methods already.

Reviewer #2: The manuscript is about exploration of what kind of biochemical PTM networks can be bistable. The author came up with a clever approach for building networks from group-up using smaller pieces (partitions) as building blocks. I personally really appreciate this topic and would like to see this manuscript published. However, it requires a lot of work on explaining all the details of the approach.

Major isssues:

The concept of partitioning isn’t really clear. It isn’t thoroughly explained in the text, nor any references were provided. I searched this topic of partitioning enzymatic networks and found nothing. I assume this is the key novelty of the paper. So this piece really requires some additional work to make sure the concept is clear to the reader.

What is the point of Figure 1? Are these building blocks of partitioning? Is this a complete set of all building blocks?

Figure 2 also leaves a lot of questions. How does Fig 2a networks leads to {{P1,S1},{P1,S2,E1,E4 },{P4,S3},{P3,S4,E3,E2}} partitioning? Same pertains to even a more complicated Figure 2b. Are partitionings unique?

I’d suggest perhaps coming up with a figure that shows in details how the network is constructed from partitions.

Is the terms “partition” the right one? Usually this means that something was divided into smaller non-overlapping parts. I think the key concept of the paper is construction of the networks from smaller blocks. This slightly non-conventional use of terminology could be confusing for the reader.

Introduction of the distance in the partition space (equation 3) would also benefit from more explanations. I would suggest using a small example and show how the equation (3) can be applied. Evaluation of bistability described in Figure 3 isn’t clear either.

I wonder about the term “controllable bistability”. The author gives the definition and refers to the citation 20. Frankly, this is the first time I see the terms controllable bistability. Is there an uncontrollable one in the chemical reaction networks? Please provide an example of uncontrollable bistability so the use of “controllable” part is justified.

Results section, step 4 of the algorithm is the key for evaluating the bistability of the constructed network. What is the “degree” of controllable bistability? Then algorithm itself is kind of brief. Evaluation for bistability of a parameter-free network is a huge problem. I don’t fully understand what compromises have been taken in step 4, so the evaluation is quite short. It feels like it was a direct simulation of ODEs to reach the equilibrium and then concentration of output species plotted against the input. If this is the case, then it is a very risky approach. How to be sure that the hyperparameter space was searched exhaustively enough?

As I mentioned before, I’d love to see this published and used by the community. However, if the paper is published as is, it will be lost because it is hard to follow. If I may suggest an example of a paper that is written like an easy to follow tutorial is “Chemical Reaction Network Theory elucidates sources of multistability in interferon signaling” by Irene Otero-Muras et al, PLOS Comp Biol.

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005454&rev=2

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

24 Jun 2021

I revised the manuscript by taking account of the reviewer’s relevant comments. I tried to revise the manuscript with all my efforts and would like to give a point-by-point response to comments below. The reviewer’s comments are written in italics.

Answers to the comments by Reviewer #1:

The paper is not intelligible the way it is currently written. It is difficult to say whether the findings merit publication if written better, because in the current form it is difficult to understand even what the main point is. The author says "proposed search algorithm worked well" in the abstract without specifying the basis for this statement. Similarly, there is a vague and general remark at the end of the abstract "seems to be applicable to the search for dynamic properties".

The paper has been reviewed in general and additions have been made where necessary.

But mostly the main text of the paper is extremely unclear. Maybe the author can summarize the goal of the paper in a short paragraph, and explain one example clearly. Unfortunately, the way it is currently written, I am afraid the paper may be of little use to anyone who is not completely well-versed with the author's methods already.

The paper has been reviewed in general and additions have been made where necessary.

Answers to the comments by Reviewer #2:

The concept of partitioning isn’t really clear. It isn’t thoroughly explained in the text, nor any references were provided. I searched this topic of partitioning enzymatic networks and found nothing. I assume this is the key novelty of the paper. So this piece really requires some additional work to make sure the concept is clear to the reader.

The paper has been reviewed in general and additions have been made where necessary.

What is the point of Figure 1? Are these building blocks of partitioning? Is this a complete set of all building blocks?

The building block of partitions is only the post-translational modification reactions shown in Eq (1), and the N post-translational modification reactions are interconnected by species identification to form an enzymatic reaction network. The family of the set of chemical species to be identified is the partition representation. Fig 1 shows an example of the case where N=1 and N=2. In other words, Fig 1 is an example of an enzymatic reaction network composed of rather than building blocks.

To make the above clear, I added an example in Fig 1 for the case of N=1, and in the text, I added and modified the partition representations corresponding to each enzymatic reaction network. I also added and corrected the explanation of the network expression in the figure, which I think was missing.

Figure 2 also leaves a lot of questions. How does Fig 2a networks leads to {{P1,S1},{P1,S2,E1,E4 },{P4,S3},{P3,S4,E3,E2}} partitioning? Same pertains to even a more complicated Figure 2b. Are partitionings unique?

I’d suggest perhaps coming up with a figure that shows in details how the network is constructed from partitions.

As in Fig 1, the set M, which is the source of the partition, is explicitly shown for Fig 2. There is a one-to-one correspondence between the partition and the enzymatic reaction network that is constructed.

Is the terms “partition” the right one? Usually this means that something was divided into smaller non-overlapping parts. I think the key concept of the paper is construction of the networks from smaller blocks. This slightly non-conventional use of terminology could be confusing for the reader.

In this paper, partition is also used in the sense of a family of subsets, where the sets are disjoint to each other and their union set is the original set. In other words, it is a division into small non-overlapping parts. This is the representation of an enzymatic reaction network.

Introduction of the distance in the partition space (equation 3) would also benefit from more explanations. I would suggest using a small example and show how the equation (3) can be applied.

Concrete examples were added to the explanation of the distance between divisions (Eq 3). In addition, diagrams of the enzymatic reaction networks corresponding to the examples were added as Fig 3.

Evaluation of bistability described in Figure 3 isn’t clear either.

The method of evaluating bistability using the function in Fig 3 (Fig 5 in the revised version) has been substantially added and revised.

I wonder about the term “controllable bistability”. The author gives the definition and refers to the citation 20. Frankly, this is the first time I see the terms controllable bistability. Is there an uncontrollable one in the chemical reaction networks? Please provide an example of uncontrollable bistability so the use of “controllable” part is justified.

In Reference 20, it is called “resettable bistability”. Thus, I will change to that term as well. Futhermore, I have added Fig 4 and a description to explain the relationship between bistability and resettability.

Results section, step 4 of the algorithm is the key for evaluating the bistability of the constructed network. What is the “degree” of controllable bistability? Then algorithm itself is kind of brief. Evaluation for bistability of a parameter-free network is a huge problem. I don’t fully understand what compromises have been taken in step 4, so the evaluation is quite short. It feels like it was a direct simulation of ODEs to reach the equilibrium and then concentration of output species plotted against the input. If this is the case, then it is a very risky approach. How to be sure that the hyperparameter space was searched exhaustively enough?

I have added an explanation to the method of evaluating controllable (resettable) bistability in step 4 of the algorithm.

Also, as you pointed out, as the stable equilibrium point (steady state) of the system, I used the value at which the time evolution of the system reached equilibrium in the direct simulation of the ODE. This is because the system is nonlinear and has a high order, so solving the algebraic equations corresponding to the steady state would eventually have to be done numerically, and determining whether the system is stable or unstable would require calculations such as finding the eigenvalues of the Jacobi matrix, which I thought would not be very efficient. I will add an organized explanation of the parameters and hyperparameters of the system in the text.

As I mentioned before, I’d love to see this published and used by the community. However, if the paper is published as is, it will be lost because it is hard to follow. If I may suggest an example of a paper that is written like an easy to follow tutorial is “Chemical Reaction Network Theory elucidates sources of multistability in interferon signaling” by Irene Otero-Muras et al, PLOS Comp Biol.

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005454&rev=2

Thank you for introducing the appropriate literature. I will refer to it.

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

15 Sep 2021

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Ivan Kryven

Academic Editor

PLOS ONE

Journal Requirements:

Additional Editor Comments (if provided):

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #2: (No Response)

Reviewer #3: (No Response)

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #2: Partly

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #2: N/A

Reviewer #3: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #2: No

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #2: No

Reviewer #3: No

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #2: While the author has addressed a lot of concerns, one of my major ones remains to be fully addressed. Specifically, I don’t understand how the evaluation function for resettable bistability (section starting on page 17, line 379) works. Either I am missing some details or this evaluation must be extremely computationally intensive.

First of all, my understanding is that the author decided to take a straightforward approach by evaluating the bistability by numerically solving the generated differential equations. Let’s say even in a fairly moderate biologically-relevant ODE system, in total there are about 10 total parameters (starting material amounts and kinetic constants). It seems like the author proposed to use 21 values per dimension. This translates into 16 trillion times solving the ODE system! Even if the number of dimensions is 5 (that limits the method to really primitive systems), it requires 4 million solutions. Let’s say one evaluation takes 1 second, this translates into about 1 month of computation time. So the straightforward numerical evaluation of bistability must be cursed by high dimensionality.

Again, for the sake of the benefit of the doubt, I may be missing something in the algorithm and/or the assumptions. I’d note that evaluation for resettable or not bistability is a hard problem. There are a variety of known approaches of how to evaluate bistability. Each of those approaches has its limitations. But no one in the field even considers evaluation of bistability for direct simulation of ODEs for the reasons outlined above.

The entire “Evaluation function for resettable bistability” needs to be improved.

1. Author needs to address the dimensionality of the problem.

2. There needs to be some discussion of kinetic constants along the species concentrations.

3. What is theta in step 2 (line 394)? How to estimate this theta?

4. Why does the author need to take quadruple of the variance? Why quadruple? What is the purpose of the entire step?

5. Since that point I’ve lost track.

Also, to make sure that the manuscript complies with "Have the authors made all data underlying the findings in their manuscript fully available?” question, I’d suggest making all the code and scripts available. I understand that the manuscript is mostly theoretical. However, if the description of the algorithms is not clear, the code may help.

Reviewer #3: This manuscript introduces a novel approach to explore dynamical properties of an enzymatic network. The novelty comes from the representation of an enzymatic network as a partition of the set instead of a regulatory matrix. This representation offers a possibility for an effective search of various dynamical properties of a network. The results demonstrate the exploration of reaction networks for the properties such as bistability and resettability.

The paper demonstrates a novel approach and important results for the exploration of bistable biological systems. However, the way this paper is written is not yet suitable for publication. “Introduction” and “Materials and Methods” sections were not comprehensible to me when reading the paper for the first time. The “Results” section clarified some confusion that appeared while reading the first part of the paper. Below are comments that might improve the readability of the paper and questions that appeared while reading it:

1. Resettability and bistability are rather central concepts in this paper. These terms appear a lot in all parts of the paper except in “Materials and Methods”. It seems like resettability and bistability drop out of that section and appear again in “Results”, making it harder for the reader to reconnect to the main story of the paper. I would recommend introducing these terms already in “Materials and Methods” next to enzymatic reaction networks which exhibit these properties. Figure 4a next to a visualization of an enzymatic reaction network which has this property would serve as a good illustrative example of what you are looking for. Then, the partition representation can be introduced for the purpose of finding networks with these properties.

2. In the subsection “Relationship between bistability and resettability” it is not clear which enzymatic network is discussed. “The relationship between bistability and resettability when the input chemical species is E1 and the output chemical species is P2 is shown in Fig 4” – it is not clear which network structure is discussed in this sentence and what are E1 and P2. Is it the network from Figures 1,2,3 or 6?

3. The algorithm descriptions in the section “Exploration of enzymatic reaction networks in the partition representation space” can be improved by either representing them as block diagrams or writing pseudo-code.

4. Section “Exploration of enzymatic networks in the partition representation space” can be framed with some intro paragraph explaining what to expect from the algorithm and some concluding paragraph summarizing the subsection. The section ends too abruptly with the steps of the algorithm. Moreover, it would create better connection with the rest of the paper if you already briefly mention how the search for resettability and bistability can be incorporated in this scheme.

5. Lines 51-53: “There, the control relationship between the cyclic reaction systems is represented by a control matrix whose elements are the cyclic reaction systems that make up the system.” – please rephrase. “Systems that make up a system” sounds confusing.

6. In the sentence in lines 75-78 you list various properties and only biochemical adaptation has a definition next to it. Please, make this consistent. You could introduce the definitions to all properties in a separate sentence, or even omit the definition and add a good reference for the interested reader.

7. Sentence in lines 79-82 has too many nested clauses, which makes it hard to grasp for the reader. Please, consider rephrasing it.

8. In the end of Introduction, term “resettability” appears rather unexpectedly. Moreover, ending the section with a definition worsens the readability of an article. Try to properly incorporate this term in the Introduction.

9. Line 93: please, rephrase the second half of the sentence and avoid “let’s think”.

10. Figures 1,2,3,6: rectangles for S and P would look better if they were white inside rather than transparent.

11. From line 124 till line 148 almost every other sentence starts with “In general”. Please, remove that or come up with another opening word.

12. Line 148: “dimer formation-separation reaction”

13. Line 163: “Fig 2 shows examples of more complex enzymatic reaction networks and its parftition representations” – typo in “partition”, and I believe it should be “reaction networks and their partition representations”.

14. Line 170: remove extra “the”.

15. Line 175-176: what do you mean by inactive and active MAPKKs? Maybe add some context for a reader who is not familiar with this system.

16. Line 185: “The advantage of the partition representation is that every partition is a representation…” – please, rephrase.

17. Line 197: “…the set M shown in Eq” -> “…the set M as shown in Eq”

18. Sentence in line 234: please expand or rephrase for better readability. The introduction of sk is a bit confusing.

19. Line 270: rk -> rk

20. Line 301: “…φ(P,ν) is a function of the partition P, with the parameter ν a pair of reaction…” – I believe you meant “… the parameter ν, which is a pair of reaction …”

21. Line 309: “and takes a maximum value among them as a value” -> “and takes a maximum among them as a value”. Line 309-310 – please, rephrase and reduce the number of words “value”.

22. Sentence in lines 317-320 is too long and might confuse the reader. If I understood it correctly, it can be split into 2 sentences: “The partition search function ψ(P0, P, σ ) described below is a recursive search function whose arguments are the initial partition P0, the current partition P, and the search depth σ. ψ(P0, P, σ ) returns the partition whose evaluation function value is greater than or equal to Ф(P0).”

23. Line 325: “is greater than of P0 ”

24. Line 330: the last step says “return P0 as a value” without any condition, which sounds confusing. Please add a condition, smth like “if all steps above are completed…” or “if none of the conditions above are satisfied…”

25. Line 336: “Resettablility is the property”. Resettable is an adjective.

26. Lines 371-372: “This is an example of a bistable that is not resettable” -> “This is an example of a bistable system (or a bistability?) that is not resettable”

27. The sentence in lines 508-513 is too long. Please, rephrase for better readability.

28. Figure 7 – why do some dots have different colors than others?

29. Figure 8 – as I understood, some arrows (like (a) and (g), and (b) and (h)) overlap significantly. Could you come up with a better representation? Maybe highlight it in the description of the figure or use dashed lines.

30. Lines 558-563 seem like an outlook. Consider first summarizing your findings and writing an outlook in the end of the Discussion. Also, line 563: I -> We, please, be consistent.

31. Line 569-570: “Good results” – please, avoid using this in a scientific paper.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

17 Oct 2021

Re: Letter to Editors

October 18, 2021

PLOS ONE

Editorial Office

Dear Editors:

Please find enclosed a revised manuscript entitled “The partition representation of enzymatic reaction networks and its application for searching bi-stable reaction systems” [PONE-D-21-10250R1] - [EMID d593f28512f8cd09], for publication in PLOS ONE.

I revised the manuscript by taking account of the reviewer’s relevant comments. I tried to revise the manuscript with all my efforts and would like to give a point-by-point response to comments below. The reviewer’s comments are written in italics.

Answers to the comments by Reviewer #2:

While the author has addressed a lot of concerns, one of my major ones remains to be fully addressed. Specifically, I don’t understand how the evaluation function for resettable bistability (section starting on page 17, line 379) works. Either I am missing some details or this evaluation must be extremely computationally intensive.

First of all, my understanding is that the author decided to take a straightforward approach by evaluating the bistability by numerically solving the generated differential equations. Let’s say even in a fairly moderate biologically-relevant ODE system, in total there are about 10 total parameters (starting material amounts and kinetic constants). It seems like the author proposed to use 21 values per dimension. This translates into 16 trillion times solving the ODE system! Even if the number of dimensions is 5 (that limits the method to really primitive systems), it requires 4 million solutions. Let’s say one evaluation takes 1 second, this translates into about 1 month of computation time. So the straightforward numerical evaluation of bistability must be cursed by high dimensionality.

Again, for the sake of the benefit of the doubt, I may be missing something in the algorithm and/or the assumptions. I’d note that evaluation for resettable or not bistability is a hard problem. There are a variety of known approaches of how to evaluate bistability. Each of those approaches has its limitations. But no one in the field even considers evaluation of bistability for direct simulation of ODEs for the reasons outlined above.

The entire “Evaluation function for resettable bistability” needs to be improved.

The method of constructing the evaluation function for resettable bistability has been modified to address each of the points pointed out below, and the overall description has also been modified to make it easier to understand.

1. Author needs to address the dimensionality of the problem.

The reviewer estimates that even if there are 5 parameters, if each value is 21 ways, we get 5^21 ≈ 4000000, which is equivalent to an exhaustive evaluation on the parameter values. I used to do this kind of exhaustive analysis, but as the reviewer pointed out, it is computationally explosive, so this is not the method used in this study. In this study, only five random combinations of parameter values are used for a single enzymatic reaction network, and the evaluation is done with those values. This value of 5 is the value of the parameter λ.

It should be pointed out, however, that the same enzymatic reaction network may be revisited, as we do not try to remove enzymatic reaction networks that have been evaluated once from the evaluation. Even in this case, the search depth σ of the two-branch search with two different partition modifiers is set to 4, so the maximum number of nodes in the binary tree is 62=2^1+2^2+2^3+2^4+2^5, and the maximum number of iterations of the main loop to change the partition is 150, so the computational cost is at most 9600=62*150 times. In addition, the actual number of parameters such as reaction rate constants and total concentration is not five, but four post-translational modification reactions, so the maximum number of reaction rate constants is 12=3*4 and the maximum number of total concentrations is 8=2*4, so the total number is a maximum of 20. This means that the cost of evaluation in the parameter space of one enzymatic reaction network is 5*9600 versus 21^20.

Why such a sparse sampling of parameter values can successfully reveal an enzymatic reaction network with the desired properties may require some explanation. Based on the exhaustive analysis I have done, one explanation is that the bistable nature of enzymatic reaction networks is dominated by the network structure and is quite robust in terms of changes in parameter values. This means that the bistable property is maintained over a quite wide range of parameter values. Alternatively, we can consider that a robust enzymatic reaction network is selectively found. However, it should be noted that it is not robust with respect to input chemical species, as we have placed the condition that it is resettable.

2. There needs to be some discussion of kinetic constants along the species concentrations.

The rate constants of the post-translational modification reactions and the total concentration of enzymes in the enzymatic reaction network are set to include approximately the values reported for the kinases that make up the MAPK cascade, which is known as a typical signaling system. This description was added in the revised text.

3. What is theta in step 2 (line 394)? How to estimate this theta?

The total concentration of each enzyme is randomly assigned as the initial value to satisfy the conservation law for each enzyme, and the process of finding 21 pairs of steady state values for the output chemical species is repeated theta times. As a result, we obtain theta steady states when the reaction rate constants of the constituent post-translational modification reactions and the total concentration of the enzyme are fixed. When all these values are the same, it means monostability, and when they are divided into two kinds of values, it means bistability. This explanation has been added to the relevant section. The value is set to 5 as appropriate, but as described in the discussion in the text, we tried other values and found no significant difference in the search results.

4. Why does the author need to take quadruple of the variance? Why quadruple? What is the purpose of the entire step?

From the theta steady state values obtained up to the previous step, the variance of those values can be used to quantitatively evaluate the degree of bistability. Since the steady state values are normalized in the range of 0 to 1, the variance takes a maximum value of 1/4 in a perfect bistable state, i.e., when half of the values are 0 and the rest are 1, and a minimum value of 0 in a monostable state, i.e., when all values are the same. So we multiply this value by 4 for the purpose of normalizing it. This explanation has been added to the corresponding step.

5. Since that point I’ve lost track.

As I answered in the beginning, I tried to revise the description to make it clearer for the other parts you pointed out.

Also, to make sure that the manuscript complies with "Have the authors made all data underlying the findings in their manuscript fully available?” question, I’d suggest making all the code and scripts available. I understand that the manuscript is mostly theoretical. However, if the description of the algorithms is not clear, the code may help.

For reference, the program code is submitted as supplementary material.

Answers to the comments by Reviewer #3:

This manuscript introduces a novel approach to explore dynamical properties of an enzymatic network. The novelty comes from the representation of an enzymatic network as a partition of the set instead of a regulatory matrix. This representation offers a possibility for an effective search of various dynamical properties of a network. The results demonstrate the exploration of reaction networks for the properties such as bistability and resettability.

The paper demonstrates a novel approach and important results for the exploration of bistable biological systems. However, the way this paper is written is not yet suitable for publication. “Introduction” and “Materials and Methods” sections were not comprehensible to me when reading the paper for the first time. The “Results” section clarified some confusion that appeared while reading the first part of the paper. Below are comments that might improve the readability of the paper and questions that appeared while reading it:

1. Resettability and bistability are rather central concepts in this paper. These terms appear a lot in all parts of the paper except in “Materials and Methods”. It seems like resettability and bistability drop out of that section and appear again in “Results”, making it harder for the reader to reconnect to the main story of the paper. I would recommend introducing these terms already in “Materials and Methods” next to enzymatic reaction networks which exhibit these properties. Figure 4a next to a visualization of an enzymatic reaction network which has this property would serve as a good illustrative example of what you are looking for. Then, the partition representation can be introduced for the purpose of finding networks with these properties.

Following your suggestion, I moved "Relationship between Bistability and Resettability" to " Materials and Methods" and changed it to "Resettable Bistability". I also referred to the enzymatic reaction network of Fig 2a as the example.

2. In the subsection “Relationship between bistability and resettability” it is not clear which enzymatic network is discussed. “The relationship between bistability and resettability when the input chemical species is E1 and the output chemical species is P2 is shown in Fig 4” ? it is not clear which network structure is discussed in this sentence and what are E1 and P2. Is it the network from Figures 1,2,3 or 6?

In "Resettable Bistability", which was moved and renamed in response to comment 1, I added an explanation referring to the enzymatic reaction network of Fig 2a as the example, and also added an explanation for E1 and P2 in that part.

3. The algorithm descriptions in the section “Exploration of enzymatic reaction networks in the partition representation space” can be improved by either representing them as block diagrams or writing pseudo-code.

I tried to break down each step of the algorithm and modify it to be at the level of pseudocode.

4. Section “Exploration of enzymatic networks in the partition representation space” can be framed with some intro paragraph explaining what to expect from the algorithm and some concluding paragraph summarizing the subsection. The section ends too abruptly with the steps of the algorithm. Moreover, it would create better connection with the rest of the paper if you already briefly mention how the search for resettability and bistability can be incorporated in this scheme.

Following your suggestion, I rearranged the description of the entire section and added the intro and final paragraphs.

5. Lines 51-53: “There, the control relationship between the cyclic reaction systems is represented by a control matrix whose elements are the cyclic reaction systems that make up the system.” ? please rephrase. “Systems that make up a system” sounds confusing.

I deleted the unnecessary part of the sentence you mentioned.

6. In the sentence in lines 75-78 you list various properties and only biochemical adaptation has a definition next to it. Please, make this consistent. You could introduce the definitions to all properties in a separate sentence, or even omit the definition and add a good reference for the interested reader.

I added explanations and references for all five properties mentioned.

7. Sentence in lines 79-82 has too many nested clauses, which makes it hard to grasp for the reader. Please, consider rephrasing it.

I broke the sentence you pointed out into two sentences to eliminate nesting, and also modified the next sentence slightly.

8. In the end of Introduction, term “resettability” appears rather unexpectedly. Moreover, ending the section with a definition worsens the readability of an article. Try to properly incorporate this term in the Introduction.

Resettability was introduced as an additional property of bistability, and an explanation was added.

9. Line 93: please, rephrase the second half of the sentence and avoid “let’s think”.

I paraphrased using "Considering".

10. Figures 1,2,3,6: rectangles for S and P would look better if they were white inside rather than transparent.

I followed your suggestion and whitened the inside of the rectangle.

11. From line 124 till line 148 almost every other sentence starts with “In general”. Please, remove that or come up with another opening word.

I’ve removed all "in general."

12. Line 148: “dimer formation-separation reaction”

I have corrected it as you suggested.

13. Line 163: “Fig 2 shows examples of more complex enzymatic reaction networks and its parftition representations” ? typo in “partition”, and I believe it should be “reaction networks and their partition representations”.

You are correct. I’ve corrected it.

14. Line 170: remove extra “the”.

I’ve removed the extra “the”.

15. Line 175-176: what do you mean by inactive and active MAPKKs? Maybe add some context for a reader who is not familiar with this system.

I followed your advice and added a brief description of the MAPK cascade.

16. Line 185: “The advantage of the partition representation is that every partition is a representation…” ? please, rephrase.

I paraphrased that sentence.

17. Line 197: “…the set M shown in Eq” -> “…the set M as shown in Eq”

I’ve fixed the sentence.

18. Sentence in line 234: please expand or rephrase for better readability. The introduction of sk is a bit confusing.

I rewrote the text, including the introduction of sk.

19. Line 270: rk -> rk

I’ve corrected the part you pointed out.

20. Line 301: “…φ(P,ν) is a function of the partition P, with the parameter ν a pair of reaction…” ? I believe you meant “… the parameter ν, which is a pair of reaction …”

You are correct. I’ve corrected the sentence.

21. Line 309: “and takes a maximum value among them as a value” -> “and takes a maximum among them as a value”. Line 309-310 ? please, rephrase and reduce the number of words “value”.

I revised the sentences and reduced the number of the word "value".

22. Sentence in lines 317-320 is too long and might confuse the reader. If I understood it correctly, it can be split into 2 sentences: “The partition search function ψ(P0, P, σ ) described below is a recursive search function whose arguments are the initial partition P0, the current partition P, and the search depth σ. ψ(P0, P, σ ) returns the partition whose evaluation function value is greater than or equal to Ф(P0).”

As you suggested, I split the sentence into two.

23. Line 325: “is greater than of P0 ”

I’ve fixed the sentence.

24. Line 330: the last step says “return P0 as a value” without any condition, which sounds confusing. Please add a condition, smth like “if all steps above are completed…” or “if none of the conditions above are satisfied…”

I have added the condition as you suggested.

25. Line 336: “Resettablility is the property”. Resettable is an adjective.

I’ve corrected the part you pointed out.

26. Lines 371-372: “This is an example of a bistable that is not resettable” -> “This is an example of a bistable system (or a bistability?) that is not resettable”

I’ve corrected the part you pointed out.

27. The sentence in lines 508-513 is too long. Please, rephrase for better readability.

I’ve corrected the part you pointed out to make it easier to read.

28. Figure 7 ? why do some dots have different colors than others?

The colors correspond to different initial values for obtaining the steady state values. However, there is no point in showing it in this figure, so we changed it to the same color.

29. Figure 8 ? as I understood, some arrows (like (a) and (g), and (b) and (h)) overlap significantly. Could you come up with a better representation? Maybe highlight it in the description of the figure or use dashed lines.

As you pointed out, some of the arrows overlap in Fig 8a. In particular, the evaluation value reaches 0.6 in all cases up to the first two steps, and almost everything else overlaps up to that point. So I changed the range of the vertical axis to 0.55 or higher, and modified it to reduce the overlap as much as possible after that. I added a note about this correction in the text and in the explanation.

30. Lines 558-563 seem like an outlook. Consider first summarizing your findings and writing an outlook in the end of the Discussion. Also, line 563: I -> We, please, be consistent.

Following your suggestion, I have revised the order to summarize the results first, and then the outlook at the end. The subject in line 563 has also been changed to "We".

31. Line 569-570: “Good results” ? please, avoid using this in a scientific paper.

I’ ve corrected the description to be more specific.

We hope the revised paper is of interest for the readers of PLOS ONE. I am looking forward to hearing from you again.

Yours sincerely,

Takashi NAKA

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

13 Jan 2022

The partition representation of enzymatic reaction networks and its application for searching bistable reaction systems

PONE-D-21-10250R2

Dear Dr. Naka,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Ivan Kryven

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #4: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #4: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #4: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #4: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #4: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #4: The work by Dr. Naka demonstrates a new computational approach to finding bistability in biochemical reaction networks. The work is technically sound and performed on an excellent professional level. The author adequately addressed comments from previous rounds of revision. I recommend this manuscript for publication at this stage.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #4: No

17 Jan 2022

PONE-D-21-10250R2

The partition representation of enzymatic reaction networks and its application for searching bi-stable reaction systems

Dear Dr. Naka:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Ivan Kryven

Academic Editor

PLOS ONE