Robustness and parameter geography in post-translational modification systems.

Biological systems are acknowledged to be robust to perturbations but a rigorous understanding of this has been elusive. In a mathematical model, perturbations often exert their effect through parameters, so sizes and shapes of parametric regions offer an integrated global estimate of robustness. Here, we explore this "parameter geography" for bistability in post-translational modification (PTM) systems. We use the previously developed "linear framework" for timescale separation to describe the steady-states of a two-site PTM system as the solutions of two polynomial equations in two variables, with eight non-dimensional parameters. Importantly, this approach allows us to accommodate enzyme mechanisms of arbitrary complexity beyond the conventional Michaelis-Menten scheme, which unrealistically forbids product rebinding. We further use the numerical algebraic geometry tools Bertini, Paramotopy, and alphaCertified to statistically assess the solutions to these equations at ∼109 parameter points in total. Subject to sampling limitations, we find no bistability when substrate amount is below a threshold relative to enzyme amounts. As substrate increases, the bistable region acquires 8-dimensional volume which increases in an apparently monotonic and sigmoidal manner towards saturation. The region remains connected but not convex, albeit with a high visibility ratio. Surprisingly, the saturating bistable region occupies a much smaller proportion of the sampling domain under mechanistic assumptions more realistic than the Michaelis-Menten scheme. We find that bistability is compromised by product rebinding and that unrealistic assumptions on enzyme mechanisms have obscured its parametric rarity. The apparent monotonic increase in volume of the bistable region remains perplexing because the region itself does not grow monotonically: parameter points can move back and forth between monostability and bistability. We suggest mathematical conjectures and questions arising from these findings. Advances in theory and software now permit insights into parameter geography to be uncovered by high-dimensional, data-centric analysis.

• We believe the reviewer has already pointed out one of our significant findings: the importance of assuming realistic enzyme mechanisms for drawing conclusions about biological robustness. As we noted in the paper, there is widespread dependence in the literature on the Michaelis-Menten enzyme mechanism. The fact that this is unrealistic because of product rebinding has been known for years, with apparently no effect. Here, we show in a compelling way why it can be so misleading. To say, as the reviewer does, that "bistability is less likely to occur", does not draw out the biological implications. Our results show that the probability of bistability reduces from 20% under the strongly irreversible Michaelis-Menten mechanism to only 1.1% for any realistic mechanism. The former figure suggests that bistability may be quite robust to change in parameter values; the latter figure does not support this conclusion at all. We believe this is a significant contribution to understanding the robustness of PTM systems. To address the reviewer's concern, we have added a paragraph to the Discussion (page 22, lines 913-9 17 ) which explains more clearly the significance of our results for biological robustness.
The reviewer described the finding above as the "major contribution" of the paper, after its "algorithmic approach". We would also draw attention to the two mathematical conjectures and question which we state in the Discussion. We believe these are also important contributions which suggest new directions of research into biochemical networks.
2. "Can we understand it mechanistically, why bistability requires that substrate and enzyme concentrations are like described in the paper?" • We were uncertain as to what precisely the reviewer is referring to here but the question seems to concern the existence of a threshold which must be exceeded for bistability to occur, as described in our first Conjecture in the Discussion section. We agree that it is interesting to ask what underlies this threshold "mechanistically" but this seems beyond the scope of the present paper. The contribution of the paper has been to reveal the existence of the threshold and to conjecture the form it should take for more general systems. We feel that we must leave it to subsequent work to explain the threshold in more detail.

Reviewer #2.
We are pleased that the reviewer found that "overall, the paper is well written", "overall this is a good balance between application, theory and computation" and recommends "it to be published".
3. "The visibility ratio sounds like it is approximating geodesics. Can numerical algebraic geometry approximate this variety in parameter space, and provide some indication of the degree of this "visibility" curve?" • We are not sure what the reviewer has in mind here. The visibility ratio is a number which measures how close the region is to being convex. There is a concept of "geodesic convexity", which could be considered for the surface of the region, with the induced Riemannian metric coming from the standard metric on the ambient parameter space. However, this is a different idea to the normal convexity of the region, for which the visibility ratio was calculated. As for the second question, we are using numerical algebraic geometry to identify points in the region and to estimate the visibility ratio. We are not sure what "visibility curve" the reviewer has in mind here, whose degree is being sought, as we have not seen this concept in the literature.
4. "The ``grammar" sounds similar to studying the graph rather than the resulting equations, which is a feature of chemical reaction network theory. It would be useful to explain the connection when presenting chemical Eqs (1). For example, this grammar sounds similar to generalized mass action kinetics." • Indeed, the grammar is formulated in the linear framework, which, as we stated (previous line 122), is "graph-based". The linear framework is related to chemical reaction network theory, as noted in the original publication (citation 69), but the corresponding graphs have different interpretations. We feel it would be too much of a digression to explain this in the text but we have reiterated the graph aspects in the relevant paragraph (lines 134-135). There is no connection to generalised mass-action kinetics and ordinary mass action is used throughout.
5. "The authors should discuss how their parameter geography compares to the work below, and whether the methods described below could provide insight to prove their conjectures:" • Thank you! We cited Conradi, Feliu et al previously (citation 58, now citation 56) but we had not taken note of papers by Joshi & Shiu and Harrington, Mehta et al. We hope to have identified these correctly from what the reviewer writes as new citations 60 and 67, respectively. We greatly appreciate these references. We have mentioned these works in the Introduction (lines 75 -78 ) and added a paragraph to the Discussion (lines 96 1 -9 66 ) , in which we review how these and other approaches may provide insight into our conjectures.
6. "For the parameter values that could not be certified, are these near the discriminant locus?" • There are different types of "discriminant locus" depending on what kind of non-generic behaviour is being discriminated. Harrington, Mehta et al discuss some of the choices (citation