Applying an inverse homeostasis perspective to simplify the design and implemention of robustly- nearly -homeostatic biological networks

A nearly-homeostatic system is able to maintain the steady state system output close to a ﬁxed set point, despite varying levels of environmental perturbation. Studying such systems enables the design of synthetic cellular systems able to function consistently under a wide range of environmental conditions. Here we present the inverse homeostasis perspective, a novel approach to studying and designing nearly-homeostatic systems. It represents a graphical approach allowing visualization of how each regulatory parameter aﬀects near-homeostatic performance while being more accessible than situation-speciﬁc, mathematically complex approaches. Given the diﬃculty of precise parameter measurement in biology, we focus on adjusting experimentally-determined steady state response curves to implement near-homeostatic systems without requiring explicit parameter measurements. One implication of the inverse homeostasis perspective is an illustration of the much stronger dependence of nearly-perfect integral con-trollers on the properties of the existing system being controlled, compared to the system-independence of perfect integral controllers.


Introduction
Near-homeostasis is the biological property of maintaining a nearly-fixed state despite different external perturbation levels. (By contrast, perfect homeostasis is the property of maintaining an exactly fixed state in the face of external perturbations, an idealization rarely realizable in practice.) Near-homeostasis is generally achieved by coupling a biological sensing mechanism to a feedback loop, so when the system experiences a new persistent perturbation level, the feedback mechanism acts to resist these changes and nearly restore the system to its original state [1,2] (Figure 1a). Nearly-homeostatic systems are abundant in biology.
In physiology, near-homeostasis is often discussed in terms of keeping global physiological parameters like internal temperature nearly constant, despite changes in environmental stimuli like external temperature that would otherwise greatly perturb those physiological parameters [3,4,5,6]. In the microbial world, a classic example is bacterial chemotaxis, where bacteria follow nutrient concentration gradients by exhibiting temporal sensory imperfect adaptation to these concentrations, responding only to changes and nearlyadapting to unchanging concentrations [7,8]. Furthermore, nearly-homeostatic systems in biology not only produce a narrow range of steady state output under wide range of perturbation levels, but they can do so even when regulatory parameters (e.g. enzyme affinity or maximum activation rate) of the system change within some particular range, a property we and others call robust near-homeostasis [9,10] (Figure 1b).
(The narrow steady state output range can change depending on the regulatory parameters of the system.) One goal in synthetic biology is to produce cellular systems that can function consistently under a range of environmental conditions, leading to an interest in the ability to design robustly nearly-homeostatic controllers to keep key system states nearly constant in the face of external perturbations. (Imagine, for example, an engineered bacterium designed to operate consistently in the human intestine, while exposed to varying extracellular conditions in terms of nutrient and ionic concentrations). While some mathematically guided robustly nearly-homeostatic designs have been proposed [11,9,12,13,14,15] and successfully implemented experimentally [16,17,15,18,19,20,21], there is no unifying theory of near-homeostasis that allows us to systematically assess the impact of each regulatory parameter on the system's ability to resist perturbations and maintain its output within a desired range. In all instances of successfully implemented robustly-nearlyhomeostatic systems [17,15,18,19], the theoretical insight of how certain regulatory parameters enhance homeostatic performance was used to enhance homeostatic performance of the experimental system. Unfortunately in these studies, the parameter insights applied were experimental-system specific; and the insights were derived using a variety of situation-specific methods that required advanced mathematical knowledge.
In addition, since many biological regulatory parameters are difficult or impossible to measure directly, if a design approach is to provide practical guidance in synthetic circuit design, it should enable us to achieve  Figure 1: Biological homeostasis. (1a) Homeostasis and near-homeostasis. A biological system generates an "output" (where we use this term to include any state of interest inside the system) in response to an applied external perturbation (hormonal signal, for example). The system is homeostatic if a shift in this perturbation from persistent level I 1 to persistent level I 2 results in the system output settling back to its original steady state(O 2 = O 1 ). That system is nearly homeostatic if a shift in this perturbation from persistent level I 1 to persistent level I 2 results in the system output settling back to near its original steady-state (O 2 ≈ O 1 ). The dashed lines illustrate that O 2 = O 1 .(1b) Robust homeostasis and robust near-homeostasis. Suppose the system with regulatory parameter set p A is homeostatic at O 2A = O 1A . The system is robustly homeostatic if for each parameter set p B in a region (more precisely, an open set) containing p A , the new system with regulatory parameter set p B is able to maintain homeostasis at O 2B = O 1B . The new homeostatic set point O 1B is not necessary the same as O 1A . A nearly-homeostatic system with regulatory parameter set p A is robustly-nearly-homeostatic if for every set of parameter values p B in a region containing p A , the new system with regulatory parameter set p B is still nearly homeostatic.
desired system behaviours without requiring complete information on the underlying regulatory parameters.
The connection between near-homeostasis and perfect homeostasis is more complex than it initially appeared. It has generally been assumed that an understanding of near-homeostasis-supporting parameter regions (ranges of system parameters that lead to the desired system behaviour) could be obtained in a straightforward manner by understanding the parameter regions supporting perfect homeostasis, then relaxing some idealized assumptions about the parameters to obtain near-homeostatic behaviour. For example, consider the case of zeroth-order degradation of the system's controller species [11,9,2], in which the degradation rate of the controller species in the idealized scenario is treated as completely independent of its own concentration by fully saturating the degradation machinery associated with the controller species. If the controller species's production rate depends only on the system output level, then the rate of change of the controller species (the sum of its production rate and degradation rate) depends only on the system output level, and importantly does not depend on the perturbation level. The system is perfectly homeostatic for for all perturbation values that allow the system to reach a steady state (where the rate of change for each system component is zero). In reality, however, it is only useful to describe degradation machinery of the controller species as being increasingly saturated but never completely saturated. In such a degradation regime, the controller species is described as being nearly but not completely independent of its own concentration and is predicted to be nearly homeostatic under the traditional assumption.
Under this traditional assumption, we have previously derived a sufficient condition for perfect homeostasis that not only unified prior published examples of perfect homeostatic supporting parameter regions under one theory, but also allowed us to systematically construct novel perfect-homeostasis supporting parameter regions [22]. However, we found that only 25% to 75% of the parameter sets supporting perfect homeostasis were successfully translated into near-homeostasis supporting ones [22], depending on the pattern of regulatory interactions (also called the network topology). Returning to the zeroth-order degradation example, if the controller species's production rate depends only on the system output level and the controller's degradation rate is nearly but not completely independent of its own concentration, then the system is not necessarily nearly-homeostatic. Our prior intuition was that steady states of a nearly homeostatic system should operate in a similar regime as steady states of a perfectly homeostatic system, so it was surprising that degradation being seemingly nearly saturated but not perfectly saturated is not sufficient to result in near-homeostasis. Since the ability to construct robust nearly-homeostatic systems does not trivially follow from the ability to construct robust perfectly-homeostatic systems and we could not pinpoint an obvious cause, we have pursued the current study on the characteristics of robustly nearly-homeostatic systems.
In this work, our first aim is to develop an unified approach to systematically determine the role of each system parameter in defining the perturbation range permitting near-homeostasis, in contrast to the situation-specific approaches employed in prior literature [17,15,18,19]. Our second aim is to create experimentally-realizable methods of implementing a near-homeostasis-supporting parameter set without being able to measure and target parameter values. We will accomplish our first aim by adopting an alternative to previous methods of translating from perfectly-homeostatic to nearly-homeostatic parameter sets. We call this alternative the inverse homeostasis perspective: rather than determining how a set of perturbations produces a narrow range of steady state system outputs, we instead (1) start with the desired range of outputs, (2) compute the range of perturbations that the system can reject without leaving the target output range, and (3) identify the key system parameters required to enlarge that range of rejected perturbations. The inverse homeostasis perspective has the advantage of being graphically-based and thus more accessible to non math-specialists. We use experimentally measurable steady state response curves to determine whether the current version of the experimental system displays the desired set of behaviours, and identify methods of shifting the response curves in ways that will enable the system to be robustly nearly-homeostatic. Using steady state response curves as a gauge of progress toward robust-homeostasis does not require measuring, estimating, or precisely targeting specific parameter values. We have also found that unlike in idealized robustly perfectly-homeostatic systems, in robustly nearly-homeostatic systems it is not generally possible to completely decouple the properties of the controller that acts upon the existing system from the properties of the existing system itself. That is, near-homeostatic performance's dependence on the existing system in presence of a near-integral controller is far greater than (perfect-)homeostatic performance's dependence on the existing system in presence of a perfect integral controller.

Problem formulation
We will first make precise the meaning of designing and implementing a robustly nearly-homeostatic system. Suppose we have an existing biological system consisting of a set of system components that regulate one another. The (regulatory) interactions between different system components are usually qualitatively described by terms such as positive, negative, or none (Figure 2a, 2c, 2e). More precisely, we suppose that the interacting components are described by a set of first order differential equations, where each equation describe how the rate of change of a system component depends on the level of both itself and the levels of other system components (Figure 2b, 2d, 2f). The mathematical characterization of this dependence consists of persistent variables or parameters representing various aspects of how that system component is regulated. These aspects can include activating effector affinity toward the current system component (K 12 , K 21 in Figure 2b, d/n in Figure 2f), maximum production rate (V 1 , V 2 in Figure 2b) and maximum degradation rate ( V 1 , V 2 in Figure 2b) of the current system component, as well as the external environment perturbation affecting the current system component (u in Figure 2b, 2d and 2f). Starting from an initial set of system component levels, the set of first order differential equations allows us to calculate the level of each system component as a function of time, and to evaluate the level of each system component once the rates of change of all system components have converged to zero (also called the steady state of the system component).
We will define one system component of interest as the system output (x 1 in Figure 2a, 2c and 2e) and one parameter of the system as the external perturbation (u in Figure 2a, 2c and 2e). We will only discuss the non-trivial case, where changing the persistent external perturbation level also changes the system output value at steady state. We then define a desirable range for the system output at steady state (called a set point range), as well as a desired range of persistent perturbation values. Our goal is to experimentally implement an additional regulatory system (called a controller) such that the existing system combined with the additional regulatory system is able to accommodate perturbations within the desired range without the steady state system output leaving the desired range. That is, whenever a persistent perturbation value is from the desired range, the steady state system output is within the desired set point range as well.

Design and implementation of robustly nearly homeostatic systems
The process of designing and implementing a robustly nearly-homeostatic system can be broken into four main steps: 1. Given the design constraints of the problem, select a feasible pattern of regulatory interactions (the network topology) to be used to implement the control over the system's output. Formulate a mathematical model using first order differential equations to capture the appropriate qualitative relationships in the system (this may include elements such as activation, repression, competitive inhibition, and so on; there are well-established representations available for all common modes of regulatory interaction).
Based on the network topology and its mathematical representation, draw on previous literature to identify parameter regions that will support perfect homeostasis for the system of interest.
2. Using the inverse homeostasis perspective, find the associated robust-near-homeostasis supporting parameter regions corresponding to the robust-homeostasis supporting parameter regions.
3. Experimentally implement an initial realization of the system, using insights from the inverse homeostasis perspective to ensure that the system has at least the capacity to be nearly homeostatic.
(a)ẋ x 1 activates regulator x 2 and x 2 represses x 1 . The first rate equation describes how the rate of change of x 1 is affected by perturbation u, regulator x 2 , and itself. The production rate of x 1 is linearly proportional to perturbation level u, and V 1 u can be seen as the perturbation dependent maximum transcription rate of x 1 . H 12 describes the fraction of promoter that is not blocked by the repressor x 2 : as x 2 increases, the fraction of active promoter goes from 1 to minimum basal fraction of b 12 . The maximum degradation rate of x 1 is V 1 , and the fraction of the maximum degradation rate is a positive saturating function of x 1 . The second rate equation describes how the rate of change of regulator x 2 is affected by the system output x 1 as well as itself. The production rate of x 2 is a positive saturating function of x 1 ; the degradation rate of x 2 is a positive saturating function of itself. All symbols except (x 1 , x 2 , u) are positive parameters. (2c and 2d) Schematic and mathematical description of the three component antithetic integral control motif. Perturbation u decreases amount of system output x 1 ; x 1 activates regulator z 2 ; regulator z 2 and regulator z 1 annihilate each other; z 1 activates system output x 1 . All symbols except (x 1 , z 1 , z 2 , u) are positive parameters. (2e and 2f) Schematic and mathematical description of the two component antithetic integral control motif that is equivalent to the three component antithetic integral control motif (2c and 2d). At steady state, rate of change of z 1 is zero, so the steady level of the intermediate regulator z 1 is only dependent on z 2 :z 1 = w/n d/n+z2 . Since we only care about how the system output respond to persistent perturbation at steady state, in which all rate equations equal to zero, we can replace every occurrence of z 1 by w/n d/n+z2 and drop the z 1 rate equation. This allows us examine how z 2 regulates system output x 1 and remove consideration of intermediate regulator z 1 from the system. The reduced antithetic integral feedback motif now resembles the two component transcriptional network (2a and 2b), with an additional linear degradation term in rate equation of regulator z 2 compared to x 2 .
Characterize the behaviour of this system using experimental measurements of steady state curves of key system components in response to their regulators. 4. The intersections of the steady-state response curves obtained in the previous step represent the steady states of those key system components in the full closed-loop system, but the location of these intersections may not fall within the desired range of output values when the system is presented with the range of perturbations that it is meant to reject. To correct this, we again employ the insights gained from an inverse homeostasis perspective to identify which qualitative adjustments to the parameters will shift the steady state response curves in the correct direction, without requiring quantitative knowledge of the actual parameter values themselves.
We will discuss each step in detail, below.
Step 1: Define a network topology, the associated mathematical model, and find a robusthomeostasis supporting parameter region Our first step is to find a controller architecture, consisting of additional system components that regulate the existing system components, as well as a range of system parameter values (also called the feasible parameter region) such that whenever persistent perturbation values are in the desired range, the steady state system output is within the desired set point range as well. More broadly, if a range of system parameter values traps the steady state system output in some narrow range that is not necessarily the desired set-point range under persistent perturbations from some much wider range, then we call this range of parameter values a near-homeostasis supporting parameter region. One can either pick a suitable controller architecture and its associated robust-homeostasis supporting parameter region from a large number of previously reported examples [23,24,11,9,2,14,12,10,25,26], or systematically derive possible controller architectures and their robust-homeostasis supporting parameter regions them by applying well-characterized principles of robust-homeostasis [22,27].
Step 2: Apply the inverse homeostasis approach to translate the robust homeostasis supporting parameter region into a robust near-homeostasis supporting region The second step translates each parameter set in the robust-homeostasis supporting parameter region into a robust-near-homeostasis supporting parameter set. A previously unsuccessful translation approach involved making ideal interaction characteristics subjectively close to ideal. Intuitively, it may appear that this translation could be accomplished by starting with idealized interaction characteristics, then create a system that is subjectively close to the idealization. Our initial attempts along these lines revealed that this approach is The degradation rate of x 2 shows little dependence on its own value. In the idealized scenario where x 2 degradation rate is completely independent of its own value, the x 2 rate equation satisfies V 2 H 21 (x 1 ) − V 2 = 0 at steady state, which results in a perturbation-independent system output steady state value x 1 . In the non-idealized scenario of pseudo-zeroth-order degradation of x 2 at steady state, a small fold increase in x 1 steady state from x L 1 to x U 1 is induced by a small fold increase in scaled production rate of x 1 (V 2 H 21 (x 1 )/ V 2 , represented as horizontal dashed lines in the top panel) or equivalently x 2 's fraction of maximum degradation rate ( θ(x 2 ), plotted as the black curve in the top panel), which in turn is induced by much larger fold increase in x 2 steady state from x L 2 to x U 2 (top and bottom panel: vertical dashed lines). The large fold increase in x 2 steady state is induced by a similar-sized decrease (bottom panel: horizontal dashed lines) in the fraction of active promoter unbound by x 2 (H 12 (x 2 ), plotted as the black curve in the bottom panel) or equivalently the perturbation adjusted scaled degradation rate of , are simply rearrangements of steady state system rate equations of x 2 and x 1 , respectively. Since there is a small fold increase in x 1 steady state, the large fold decrease in V 1 θ(x 1 )/(V 1 u) must be induced by an even larger fold increase in perturbation u from u L to u U . The perturbation range [u L , u U ] is called the theoretical perturbation range permitting near-homeostasis for set point range On the other hand, suppose the degradation rate of x 2 at steady state depends more on its own concentration. The same fold increase in x 1 is induced by the same fold increase in the scaled production rate of  In the idealized scenario where first order dilution is completely absent (d = 0), the regulator z 2 rate equation satisfies θx − w = 0, resulting a perturbation u-independent system output x 1 value at steady state. Suppose now d > 0. In the top panel, the black curve represents scaled degradation rate of z 2 , 1 θ w z2 d/n+z2 + dz 2 , and is equivalent to value of the system output x 1 at steady state. The scaled z 2 degradation rate is further decomposed into the scaled saturable degradation term w θ z2 d/n+z2 (gray curve) that serves as the lower bound of the scaled z 2 degradation rate, and upper limit of the scaled z 2 degradation rate w θ + d θ z 2 (straight gray line). In the bottom panel, the black curve represents scaled production rate of x 1 ( kw d+nz2 /( w θ ) = kθ d+nz2 ), and is equivalent to scaled degradation rate of x 1 ((y + d + w)x 1 /( w θ )) at steady state. Just like in 3a, one can follow the dashed lines from the top panel to the bottom panel to see how a small fold increase in steady state x 1 from x L 1 below w/θ to x U 1 above w/θ (top panel: horizontal dashed lines) can be induced by a much larger fold increase in steady state z 2 from z L 2 to z U 2 (top and bottom panel: vertical dashed lines). The large fold increase in steady state z 2 is in turn induced by a large fold decrease in scaled degradation of , which is ultimately induced by a large fold decrease in perturbation u from u U to u L . [u L , u U ] is called the theoretical perturbation range permitting near-homeostasis for the set point range On the other hand, if both steady states x L 1 and x U 1 are far greater than w/θ, then degradation of z 2 at steady state is dominated by first order degradation (dz 2 ). An increase in steady state x 1 (top panel: horizontal dashed lines) is only induced by a similarly small fold increase in steady state z 2 (top and bottom panel: vertical dashed lines). The small increase in steady state z 2 is equivalent to a small fold decrease in the scaled degradation of x 1 (bottom panel: horizontal dashed lines), which is ultimately induced by a similarly small fold decrease in perturbation u. not universally applicable: 25-75% of parameter sets with subjectively "near-ideal" interaction characteristics failed to show robust near-homeostasis in practice [22].
We have now developed a approach that we call the inverse homeostasis perspective, which offers a novel method of analyzing and designing robust nearly-homeostatic systems. One common way to study robustlynear-homeostasis supporting parameter spaces is to ask how a large fold change in persistent perturbation induces a much smaller fold change in the steady state system output. On the other hand, the inverse homeostasis perspective asks what properties a system must have such that a small fold change in steady state system output is induced by a much larger fold change in persistent perturbation. We start from two system output values that have the same fold difference as the desired set point range and compute the perturbations that would have resulted in each of the steady state system output values. We call those two perturbation values the theoretical perturbation range permitting near-homeostasis that is associated with the starting system output range. This computation process allows us to determine the exact role of each parameter in enlarging the theoretical perturbation range permitting near-homeostasis in terms of fold difference, placement of the exact theoretical perturbation range permitting near-homeostasis, as well as placement of the exact system output range that generated the perturbation range permitting nearhomeostasis.
Specifically, we break down each rate equation of key system components at steady state into multiple rationally determined parts, then plot both the relationship between these parts at steady state and how steady state values of key system components affect the steady state values of these parts in a single graph called an inverse homeostasis plot. While there is no fixed set of rules that determine what the rationally determined parts should be, it is often helpful to break the steady state equation of each key system component into the production rate and degradation rate. The production rate equation is plotted against one vertical axis and the degradation rate equation is plotted against the opposite but equivalent vertical axis ( Figure   3a, 3b, 4a, 4b). A set of useful inverse homeostasis plots explains how under suitable conditions a small fold difference in system output at steady state is induced by a big fold difference in persistent perturbation levels (Figure 3a, 3b, 4a, 4b).
In addition, once we have set up the inverse homeostasis plots, we can change the value of any chosen parameter, observe the corresponding shifts in the relevant rate curves, and finally observe the change in the fold difference of the theoretical perturbation range permitting near-homeostasis. In Tables 2a and 3a, we summarize the role of each parameter in changing the theoretical perturbation range permitting nearhomeostasis. In contrast to using pure numerical simulations to determine the effect of a parameter change in the context of specific values of all other parameters, the current approach is able to analytically capture the effect of that parameter change independent of specific values of other parameters.
It is important to note that we are no longer mimicking normal experimental procedures of trying to find the steady state system output based on a persistent perturbation value. Thus, the perturbation values generated from starting steady state system output values are not necessarily experimentally realizable.
Nevertheless, we know how to change parameters of the system such that theoretical perturbation range permitting near-homeostasis matches the desired range. Completion of this step allows us to find a controller architecture and a parameter region that enables the entire system to accommodate perturbations within the desired range without the steady state system output leaving the desired range, at least theoretically.
Step 3: Initial experimental implementation to ensure the system display basic properties of robust-near-homeostasis The third and fourth steps will allow us to experimentally implement one parameter set allowing the system to accommodate perturbations within the desired range without the steady state system output leaving the desired set-point range, while not requiring the ability the directly measure and therefore target parameter values. The third step is to demonstrate that the experimental system is capable of rejecting a large range of perturbations at all. This means implementing a set of parameters such that the fold difference in the calculated uncoupled theoretical perturbation range permitting near-homeostasis is larger than desired. One straightforward way to verify that the fold difference value is larger than desired is to determine the precise values of all model parameters from experimental data and calculate the theoretical perturbation range.
Despite its apparent simplicity, this approach presents substantial challenges: accurately estimating model parameters is difficult or impossible for many experimental systems. A second approach to verifying that the fold difference value is larger than desired is to experimentally implement a set of parameters such that the system is able to accommodate the desired perturbation without the system output leaving the desired setpoint range in the uncoupled manner ( Figure 5, 6). If the system is nearly homeostatic in the aforementioned uncoupled manner, then the theoretical uncoupled perturbation range permitting near-homeostasis in terms of fold difference is larger than desired.
From an undesirable steady state response curve, insights from the inverse homeostasis plots provide a list of parameters and directions of parameter change for the steady state response curve to reach a desired status (Table 2a,    : steady state x U 2 must be substantially greater than K2 in order for a small increase in system output steady state x1 from x L 1 to x U 1 to be induced by a large change in regulator steady state x2 from x L 1 to x U 2 ( Figure  3a top panel).
V2, V2: the scaled production rate of the regulator x2 at the system output x1, V2H21(x1)/ V2, must be sufficiently close to 1 in order for steady state x U 2 to be substantially greater than K2.
b21, K21: The fold increase in x1 is not equivalent to the fold increase in the scaled production rate of x2,V2H21(x1)/ V2 (horizontal dashed lines in the top panel of Figure 3a). If K21 decreases or b21 increases, the same small fold increase in x1 steady state would be induced by an decreased fold increase in the scaled production rate of x2, V2H21(x1)/ V2. As an extreme example, if b21 = 1 or K21 = 0, then there is 0 fold increase in the scaled production rate of x2 as x1 increases from x L 1 to x U 1 . Secondary facilitators of nearhomeostasis We assume that an small fold increase in system output steady state x1 from x L 1 to x U 1 is induced by a large fold increase in regulator steady state x2 from x L 2 to x U 2 in the manner described in the previous row.
b12: as b12 increases, the scaled production rate curve (H12, the black curve in the bottom panel of Figure 3a) shifts upward but is pinned at the upper left corner. The vertical dashed lines (x L 2 and x U 2 ) in the same graph panel does not change, but the horizontal dashed lines will both shift upward and also become closer together. Therefore, the theoretical perturbation range permitting near-homeostasis decreases. As an extreme example, as b12 converges to 1, the horizontal dashed lines will both converge to 1.

K12:
If b12 > 0, as K12 converge to zero, the horizontal dashed lines in Figure 3a will converge to b12 and thus the perturbation range permitting near-homeostasis will converge to 0. As K12 increases past some optimal point dependent on both b12 and x L 2 , the theoretical perturbation range permitting near-homeostasis decreases.
V1, K1, V1: these parameters do not have a dramatic effect on theoretical perturbation range permitting nearhomeostasis in terms of fold difference, but rather change the exact placement of the theoretical perturbation range permitting near-homeostasis to cover the desired range.

(a)
Parameter Biological description Experimental implication of parameter tuning

b12, b21
The basal production of the system output x1 and regulator x2 as a fraction of the maximal production rate, respectively.
The parameter value can be modulated by the promoter sequence.
As b12 increases, the dashed and dashed-dotted system output response curves in Figure 5 become more vertical, making it harder to draw the gray box overlapping both dashed and dashed curves of the same color.
In addition,increased basal expression b12 shifts the x1 steady state curve in response to x2 to higher x1 values at large regulator x2 levels. An analogous statement can be made for b21.

V1, V2
The maximum production rate of the system output x1 and regulator x2, respectively.
The parameter value can be modulated by copy number, the promoter sequence, and the ribosome binding site.
As V2 increases up to a certain point, the same small increase in x1 induces a larger fold increase in x2 at steady state. Further increased V2 prevents x2 from reaching steady state in response to x1.
In addition, increased maximum transcription V1 shifts the x1 steady state curve in response to x1 and u to higher x1 values regardless of x2 values. An analogous statement can be made for V2.

K12, K21
The amount of effector x2 and x1 needed for the production rate of x1 and x2 to be at half maximum, respectively.
The parameter value can be modulated by changing how tightly the effector binds to the promoter.
Increased K12 shifts the x1 steady state curve in response to x2 and u toward larger regulator x2 values. An analogous statement can be made for K21.

V1, V2
The maximal degradation rate of the system output x1 and regulator x2, respectively.
The parameter value can be modulated by adding protein degradation tags and changing concentration of the degradation machinery.
As V2 decreases to a certain point, the same small increase in x1 induces a larger fold increase in x2 at steady state. As V2 decreases further, x2in response to x1 will no longer reach a steady state.
Increased V1 shifts the x1 steady state response curve to smaller system output x1 values, regardless of x1 concentration. An analogous statement can be made for V2.

K1, K2
The amount of x1 and x2 for their respective degradation rate to reach half maximum.   ) saturates at a higher value. The same fold increase in the system output x1 from below w/θ to above w/θ, represented as horizontal dashed lines, becomes a larger absolute increase in x1. Since the slope of the upper gray line does not change, a larger absolute difference in the horizontal dashed lines translates to increased fold difference between the vertical dashed lines (z L 2 and z U 2 ). As w stays the same and d decreases, the upper gray line in the same panel has decreased slope and the lower gray curve reaches saturation at a smaller z2. The same fold increase in x1 from below w/θ to above w/θ is induced by large fold increase in z2. In summary, as w/d increases, the same fold change in x1 steady state (output) from below w θ to above w θ is induced by an increased fold change in z2 steady state (regulator). d, n: as the d/n ratio decreases, the same small change in x1 steady state from below w θ to above w θ is induced by increased fold change in z2 steady state. Graphically, decreased d/n makes the scaled z2 degradation rate curve w We assume that an small fold increase in x1 steady state is induced by large fold increase in z2 steady state in the manner described in the previous row.
θ, k, y: A fold change in z2 steady state is induced by same fold change in adjusted perturbation range y + d + u permitting near homeostasis, regardless of changes in θ, k, y. Changing θ, k, and y does allow the theoretical perturbation range u permitting near-homeostasis to match the desired range (Figure 4a: bottom panel).
Changing θ without changing w in the same proportion changes the optimal x1 (output) steady states that can be induced by largest fold change in z2 steady state (Figure 4a: top panel).

(a)
Parameter Biological description Experimental implication of parameter tuning k, θ First order rate constant for production of x1and z2, respectively.
The parameter value can be modulated by changing strength of activation by the affector.
As an example, increased k shifts the steady state curve of x1 in response to z2 toward larger x1 values. w The production rate of z1. Increased w increases the sharpest slope of the steady state z2 response curve as a function of x1.
Increased w also shifts the steady state curve x1 in response to z2 and u to higher system output x1 values Increased w also shifts the steady state curve of z2 in response to x1 toward smaller z2 values.
y Degradation rate constant of x1.
The parameter value can be modulated by adding a degradation tag to x1.
Increased y shifts the x1 steady state curve in response to z2 and u toward smaller x1 values.
d First order dilution rate of all system components.
The parameter can be modulated by changing growth rate of the the cell.
Decreased d increases the sharpest slope of the steady state curve of z2 in response to x1 Decreased d has limited ability to shift the steady state curve of z2in response to x1 toward smaller z2 values.
n The annihilation rate constant between z1 and z2.
Increased n increases the sharpest slope of the steady state curve of z2 in response to x1.
Increased n has limited ability to shift the steady state curve of z2 in response to x1 toward smaller z2 values. Figure 5: Simulated experimental implementation of a robustly-nearly-homeostatic transcriptional network in Figure 2a. We want our system to have steady state system output x 1 between x L 1 and x U 1 while the perturbation is between u L and u U . Starting from an initial set of parameters involved in the control of regulator x 2 (panel a: solid gray curve), we tune appropriate parameters involved in the control of regulator x 2 (b 21 , K 21 , V 21 , V 2 , K 2 ) according to Table 2a and 2b such that x 2 steady state response curve increases rapidly from x L 2 to x U 2 as x 1 increases from x L 1 to x U 1 (panel a: solid black curve). We define the set of points on the x 2 steady state response curve from (x L 1 , x L 2 ) to (x U 1 , x U 2 ) (panel a: the black box) as the permissible steady states of final nearly-homeostatic system, denoted by X P . We then move on to finding suitable parameter involved in the control of the system output x 1 . We start by searching for a satisfactory set of parameters involved in regulation of the system output x 1 (b 12 according to Table 2a, . In most if not all cases, the aforementioned condition is equivalent to being able to draw a box of the same size as the black box (panel b: the gray box) on the logarithmic portion of the graph (x 1 ≥ 3.4/α x1 and x 2 ≥ 3.4/α x2 ) that overlaps with both gray x 1 steady state response curves. The current set of parameters for the entire system must have larger than desired uncoupled theoretical perturbation range permitting near-homeostasis, completing step 3 of the design process. Lastly, we find suitable K 12 , V 1 , V 1 , K 1 (Table 2a, . The current set of system parameters allow the system to accommodate the desired perturbation range without the steady state system output leaving the desired set point range, completing step 4 of the design process. (c) Figure 6: Simulated experimental implementation of a robustly-nearly-homeostatic antithetic integral feedback motif in Figure 2e. We want our system to have steady state system output x 1 between x L 1 and x U 1 while the perturbation is between u L and u U . Starting from an initial set of parameters involved in the control of regulator z 2 (panel a: solid gray curve), we search for a satisfactory set of parameters involved in the control of regulator z 2 (d, n, w, θ according to Table 3a, 3b) such that z 2 steady state response curve increases from z L 2 to z U 2 rapidly as x 1 increases from x L 1 to x U 1 (panel a: solid black curve). We define the set of points on the z 2 steady state response curve from (x L 1 , z L 2 ) to (x U 1 , z U 2 ) (panel a: the black box) as the permissible steady states of final nearly-homeostatic system, denoted by X P . We then move on to finding suitable parameter involved in the control of the system output x 1 . We start by verifying that there exists a point ( x L 1 , z L 2 ) on the x 1 steady state response curve under minimal perturbation u L (panel b: gray dashed −−line) and a point x U 1 , z U 2 on the x 1 steady state response curve under maximal perturbation u U (panel b: gray dashed-dotted − · − line) that preserves ratio between two points in the permissible set of steady states X P (i.e. there exists points . In most if not all cases, the aforementioned condition is equivalent to being able to draw a box of the same size as the black box (panel b: the gray box) on the logarithmic portion of the graph (x 1 ≥ 3.4/α x1 and z 2 ≥ 3.4/α z2 ) that overlaps with both gray x 1 steady state response curves. The current set of system parameters must result in larger than desired uncoupled theoretical perturbation range permitting near-homeostasis, completing step 3 of the design process. Lastly, we adjust parameters associated with x 1 regulation (k, y according to Table 3a, 3b) such that such that x 1 response curves under the minimal perturbation u L (panel c: black dashed −− curve) and the maximal perturbation u U (panel c: black dashed-dotted − · − curve) intersect the z 2 response curves within the black box bounded by . The current set of system parameters allow the system to accommodate the desired perturbation range without the steady state system output leaving the desired set point range, completing step 4 of the design process. another describing the steady state level of a key regulatory component in response to different persistent levels of the system output ( Figure 5, 6: solid curves).
Note that the plots employ the log modulus scale, which is convenient for plotting experimental response curves in a manner that is consistent with the inverse homeostasis perspective and for estimating certain parameters. The log modulus scale plotted value = sign(measured value) log 10 (α|measured value| + 1), α > 0 transforms both positive and negative measured values to a hybrid logarithmic scale. Every two sufficiently large measured values dependent on α (≥ 3.4/α in absolute value) that are ten fold apart are plotted as two points that are nearly a fixed distance apart (0.9 to 1 unit apart) on the graph. The value of α is set such that measurements below the minimum absolute value of 3. Step 4: Adjust system parameters to achieve desired output range in the face of the required range of input perturbations.
The fourth and final step is to change remaining parameters of the system such that the system is able to accommodate the desired perturbation range while the steady state system output stays within the desired set-point range in the coupled manner. The ranges of those system components generates a bounded and possibly multi-dimensional box ( Figure 5, 6, box with solid black boundaries). Our task is to change parameters of the system to ensure that perturbation dependent steady state response curves in the presence of both lower and upper limits of perturbation intersect with all other steady state response curves within that bounded box ( Figure 5, 6: black dashed system output steady state response curves), using the summarized roles of various parameters in homeostatic performance (Table 2a, 2b, 3a, 3b).
Key properties of the algorithm for design and implementation of robustly nearly homeostatic systems Unlike systems with perfect integral control where all perturbations that permit a system steady state enable perfect homeostasis, the dependence of the fold perturbation range permitting near homeostasis on the existing system is far greater than the dependence of the fold perturbation range permitting perfect homeostasis on the existing system. Any part of the entire system can act to reduce the fold difference in perturbation range permitting near-homeostasis that is enlarged by other parts of the entire system. Suppose we have a system where an additional regulator attempts to ensure near-homeostasis of an existing negative feedback network consisting of the system output and a native regulator, by sensing the system output levels and acting only on the native regulator (Figure 7a). If the additional regulator is a perfect integral controller, then the perturbation range permitting homeostasis is independent of the existing negative feedback system ( Figure 7b). However, if the additional regulator is a near -integral controller, then the perturbation range permitting near-homeostasis is highly dependent on regime of activation of the native regulator by the additional regulator (Figure 7c). In this case, a small change in steady state system output is induced by a bigger change in steady state of the additional regulator; and the bigger change in steady state of the additional regulator is induced by a small change in level of the native regulator; and that small change in level of the native regulator at steady state is induced by a similarly small change in perturbation value.
Secondly, the design part of the algorithm is capable of seamlessly harnessing multiple poorly-performing feedback actions to cooperatively enhance the theoretical range permitting near-homeostasis in terms of fold difference. The inverse homeostasis perspective asks how a small change in steady state system output is induced by a much larger change in persistent perturbation. A small change in steady state system output can be induced by a bigger change in steady state of the regulator, and that bigger change in steady state of the regulator can be induced by an even larger change in persistent perturbation (Figure 8a, 8b, 8c,   8d). The fold difference in the theoretical perturbation range permitting near-homeostasis contributed by each feedback action cooperatively enlarges the fold difference in theoretical perturbation range permitting near-homeostasis. In our example, if basal expressions (b 21 , b 12 ) are zero, then one feedback action with a cooperativity coefficient of 9 and another feedback action with a cooperativity coefficient of 1 is equivalent  Figure 7: Near-integral feedback on an existing negative feedback network. (7a, 7b) Schematic and mathematical description. The existing system contains the system output x 1 that activates native regulator x A , who then in turn repress system output x 1 . The system output x 1 is perturbed positively by u. The additional regulator x B is activated by the output x 1 ; x B then further activates native regulator x A . The rationale of various parts of the rate equation is stated in Figure 2b. It is worth noting that the production rate of native regulator x A is maximum production rate multiplied by the fraction of promoters that is bound by both activators x 1 and x B . If x B is a perfect integral controller whose rate equation isẋ B = V B H B1 (x 1 )− V B , the perturbation range permitting homeostasis is dependent only on whether both x 1 and x A rates of change can become zero. (7c) Analysis of the system in Figure 7a from the inverse homeostasis perspective. In the top left panel, the black curve represents fraction of maximal degradation rate at different x B levels, θ B (x B ). Suppose degradation of regulator x B at steady state operates in near-integral-control regime: x B K B . A small fold increase in x 1 steady state from x L 1 to x U 1 is induced by a small fold increase in scaled production rate of x 1 (V B H B1 (x 1 )/ V B , represented as horizontal dashed lines intersecting the black curve in the top panel) and equivalently a small fold increase in x B fractional degradation rate θ B (x B ), which in turn is induced by a much larger fold increase in x B steady state (top left and bottom left panel: vertical dashed lines). This large fold increase in x B steady state is induced a smaller fold increase in scaled , for all except sufficiently large K AB . The smaller fold increase in scaled x A production rate is induced by the same small fold change in x A (bottom left panel: horizontal dashed lines). The small fold increase in x A is induced by a similar-sized decrease (bottom right panel: horizontal dashed lines) in fraction of promoter unbound by repressor x A (H 1A (x A ), plotted as the black curve in the bottom panel), which is equivalent to a similar-sized decrease in perturbation adjusted scaled degradation rate of x 1 ( V 1 θ(x 1 )/(V 1 u)). Thus, the near-homeostatic performance of the system depends on more than requirements of the perfectly-homeostatic system: K AB is sufficiently larger than both vertical dashed lines intersecting x B , in addition to x 1 and x A being able to reach steady state.
(a) (d) Figure 8: Coupling of multiple near-homeostasis facilitating feedback actions. (8a, 8b) Schematic and mathematical descriptions. The system output x 1 activates regulator x 2 ; the regulator x 2 represses the system output x 1 . The perturbation u activates the system output x 1 . (8c) Mechanism of robust-nearhomeostasis from the inverse homeostasis perspective. As the cooperativity coefficient n 21 increases, the fraction of x 2 promoter activated by x 1 show increased threshold response to x 1 level around K 21 (left panel: darkening of the gray curves represents H 21 (x 1 ) at increased n 21 ). As n 21 increases, a small increase in x 1 steady state from x L 1 smaller than K 21 to x U 1 larger than K 21 (left panel: vertical dashed lines) is induced by a higher fold increase in fraction of promoters activated by x 1 (H 21 (x 1 )), which is equivalent to a higher fold increase in scaled degradation of x 2 , V 2 θ(x 2 )/V 2 (left panel: horizontal dashed lines). Since the slope of the scaled degradation curve of x 2 decreases from V 2 /(V 2 K 2 ), the higher fold increase in scaled degradation of x 2 is induced by at least as much fold increase in x 2 steady state from x L 2 to x U 2 . Since increased cooperativity coefficient n 12 leads the fraction of x 1 promoters unbound by repressor x 2 to show stronger threshold dependent repressor binding around x 2 = K 12 (right panel: darkening of the gray curves represents H 12 (x 1 ) at increased n 12 ), the large fold increase in x 2 steady state is induced by an even larger fold decrease in fraction of x 1 promoters unbound by x 2 , H 12 (x 2 ) (right panel: horizontal dashed lines), as long as K 12 is between x L 2 and x U 2 . The large fold decrease in H 12 (x 2 ) is equivalent to the same large fold decrease perturbation adjusted scaled degradation rate of x 1 , V 1 θ(x 1 )/(V 1 u) (right panel: vertical dashed lines), which is induced by similarly large fold increase in perturbation u. The threshold dependent activation of x 2 by x 1 and the threshold dependent repression of x 1 by x 2 enlarge perturbation range permitting nearhomeostasis in a cooperative manner. (8d) In the absence of coupling of two threshold dependent feedback actions, like when the activation of x 2 by x 1 no longer operates in a threshold dependent manner, a small fold difference in x 1 (left panel: vertical dashed lines) is induced by a much smaller fold difference in x 2 (left panel: horizontal dashed lines). The much smaller fold difference in x 2 is ultimately induced by a much smaller fold difference in perturbation u compared to Figure 8c. In fact, if basal expressions (b 21 , b 12 ) are zero, then a single highly cooperative feedback action with a cooperativity coefficient of 9 is equivalent to two feedback actions both with cooperativity coefficient of 3.
to two feedback actions both with cooperativity coefficient of 3.

Discussion
Our previously published strong cofactor condition [22] both characterizes robust-homeostasis-supporting parameter regions and describes how such parameter regions can be constructed de novo. However, this condition does not describe how to translate from a robust-homeostasis supporting parameter set involving idealized interactions into a robust near-homeostasis supporting parameter set that is more experimentally realistic. This inverse homeostasis perspective in this study enables us to translate a robust homeostasis supporting parameter region into a robust-near-homeostasis supporting parameter region. The inverse homeostasis perspective facilitate the translation process by looking at how a small fold change in steady state system output is induced by much larger fold change in persistent perturbation and therefore giving us insight into the role of each parameter in the fold perturbation range permitting near-homeostasis. There are several additional important properties of the inverse homeostasis perspective. First, unlike the situation specific methods that require advanced mathematical knowledge, the inverse homeostasis perspective is systematic, applicable to all scenarios, primarily graphical and can be applied by non math-specialists. Second, the inverse homeostasis perspective explains why the existing system plays a far greater role in determining the near-homeostatic performance than in determining perfect-homeostatic performance. A desirable consequence of the increased dependence of near-homeostatic performance on every part of the system is that multiple poorly performing feedback actions can cooperatively enhance near-homeostatic performance.
Finally, to implement a robust-near-homeostasis supporting parameter set whose perturbation and steady state output range permitting near-homeostasis are within the desired range, we use experimental steady state response curves of key system components to give us an idea of how close we are to implementing such a parameter set and the inverse homeostatic perspective to identify parameter tuning directions to get us there. Since we only rely on parameter tuning directions, there is no need to measure any parameter values.
Although our design strategy does not require accurate parameter estimation and uses routinely applied experimental steady state response curves to implement robustly-nearly-homeostatic systems, two seemingly true assumptions of our design strategy should be experimentally verified in order to demonstrate feasibility of the design strategy. First, the overlap of experimental steady state response curve of individual system components must predict the steady state of the entire system. Second, change in experimental steady state response curve of individual system components must translate to change in perturbation dependent steady state response of the entire system (i.e. the homeostatic performance of the entire system). While the two aforementioned statements are trivially true from a modelling point of view, the two statements may not be perfectly correct in an experimental setting. For instance, one natural way to experimentally implement our two component negative transcriptional feedback system (Figure 2a) is to implement the rate equation associated with each system component as a separate transcriptional unit consisted of a regulator responsive promoter and a downstream coding sequence corresponding to that system component. In order to determine steady state response curve of a system component, the cell would contain only the relevant transcriptional unit and an auxilliary transcriptional unit that generates persistent regulator levels for that relevant transcriptional unit. The presence of the auxilliary transcriptional unit may change apparent regulatory parameters of the system component compared to the full negative feedback system. Overlap of experimental response curve of the system output with experimental response curve of the regulator might only roughly predict location of the steady state of the negative feedback system.

Author contributions
Z.T. conceived and computationally validated the inverse homeostasis perspective. Z.T. and D.R.M conceived and computationally validated the proposed approach to experimentally implement robustly-nearlyhomeostatic systems without measuring parameter values. Z.T. and D.R.M wrote the paper.