Homeostasis in networks with multiple inputs

Homeostasis, also known as adaptation, refers to the ability of a system to counteract persistent external disturbances and tightly control the output of a key observable. Existing studies on homeostasis in network dynamics have mainly focused on ‘perfect adaptation’ in deterministic single-input single-output networks where the disturbances are scalar and affect the network dynamics via a pre-specified input node. In this paper we provide a full classification of all possible network topologies capable of generating infinitesimal homeostasis in arbitrarily large and complex multiple inputs networks. Working in the framework of ‘infinitesimal homeostasis’ allows us to make no assumption about how the components are interconnected and the functional form of the associated differential equations, apart from being compatible with the network architecture. Remarkably, we show that there are just three distinct ‘mechanisms’ that generate infinitesimal homeostasis. Each of these three mechanisms generates a rich class of well-defined network topologies—called homeostasis subnetworks. More importantly, we show that these classes of homeostasis subnetworks provides a topological basis for the classification of ‘homeostasis types’: the full set of all possible multiple inputs networks can be uniquely decomposed into these special homeostasis subnetworks. We illustrate our results with some simple abstract examples and a biologically realistic model for the co-regulation of calcium (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{Ca}$$\end{document}Ca) and phosphate (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{PO}_4$$\end{document}PO4) in the rat. Furthermore, we identify a new phenomenon that occurs in the multiple input setting, that we call homeostasis mode interaction, in analogy with the well-known characteristic of multiparameter bifurcation theory.


Introduction
A homeostatic process is characterized by the following property: approximately zero steady-state error to external disturbance, which means that an observable of interest is tightly controlled.Homeostasis is biologically important because it protects organisms against changes induced by the environment.A familiar example is thermoregulation, where the body temperature of an organism remains roughly constant despite variations in its environment (Morrison 1946).Another example is a biochemical reaction network, where the equilibrium concentration of some important molecule might not change much while the concentration of another reactant changes (Reed et al. 2017).Further examples include regulation of cell number and size (Lloyd 2013), sleep control (Wyatt et al. 1999), and expression level regulation of housekeeping genes (Antoneli et al. 2018).
Homeostasis can be mathematically defined as follows (see Sect. 2.1).Consider a dynamical system depending on an external parameter I which varies over an open interval ]I 1 , I 2 [ of external stimuli.Suppose there is a family of equilibrium points X (I) and an observable φ such that the input-output function z(I) = φ(X (I)) is well-defined on ]I 1 , I 2 [.In this situation, we say that the system exhibits homeostasis if, under variation of the external parameter I, the input-output function z(I) remains 'approximately constant' over the interval of external stimuli.
There are two formulations of 'approximately constant' often considered by researchers.The first, more stringent, called perfect homeostasis, is widely studied in control engineering and synthetic biology under the name 'perfect adaptation' [cf.Mello and Tu (2003); Ma et al. (2009); Ang and McMillen (2013); Araujo and Liotta (2018); Khammash (2021); Frei and Khammash (2021)].Perfect homeostasis is defined as the ability of a system to reset to its pre-stimulated output level, called the set point, after responding to arbitrary external stimuli.It is obvious that this condition is equivalent to the requirement that the input-output function is identically constant.
The second, more general, called near-perfect homeostasis, requires that the inputoutput function stays within a 'narrow' range under variation of external stimuli over a bounded interval.Hence, a typical 'plot' of the input-output function has a bounded region of homeostasis where it is approximately constant, called the plateau, flanked by regions of escape from homeostasis, where it varies monotonically.See plots of inputoutput functions fitting data sets sampled from real biological systems in Morrison (1946), Golubitsky and Stewart (2017b), Nijhout et al. (2018).
The notion of near-perfect homeostasis has appeared in the literature under the names near-perfect adaptation (Mello and Tu 2003;Ang and McMillen 2013;Ferrell 2016) and imperfect adaptation (Bhattacharya et al. 2022(Bhattacharya et al. , 2021)).A refinement of the notion of near-perfect homeostasis, called infinitesimal homeostasis has been proposed by Golubitsky and Stewart (Golubitsky and Stewart 2017b).Since then, aspects of this new concept have been explored in several publications (Reed et al. 2017;Golubitsky and Stewart 2018;Duncan et al. 2018;Duncan and Golubitsky 2019;Golubitsky and Wang 2020;Wang et al. 2021;Madeira and Antoneli 2022).
In this paper, we shall study near-perfect adaptation from the point of view of infinitesimal homeostasis theory.Other contributions on near-perfect homeostasis include (Bhattacharya et al. 2023;Blanchini et al. 2022;Gross et al. 2019) and references therein.These three groups propose distinct approaches to near-perfect homeostasis, under various assumptions on the functional form of the dynamics.Moreover, since all these approaches are distinct from the one presented here, it would be a too large detour for us to try to explain and compare all these ideas.Nevertheless, we believe that a review unifying all these ideas would an extremely valuable achievement that could bring closer the several groups working in this subject under different viewpoints.
We say a system with input-output function z(I) exhibits infinitesimal homeostasis if dz dI (I 0 ) = 0 for some input value I 0 ∈]I 1 , I 2 [.The vanishing of the derivative of z at I 0 implies that I 0 is a critical point of z.Moreover, the second order derivative of z with respect to I 0 can be used to give a quantitative estimate on the size of the interval ]I 1 , I 2 [ where z(I) stays within z(I 0 ) ± δ, for a given δ > 0 (see Golubitsky and Stewart (2022) for details).As we shall see in a moment, there are some additional advantages in adopting this point of view, besides providing a plausible notion of near-perfect homeostasis.
Nijhout, Reed, and Best (Nijhout et al. 2004(Nijhout et al. , 2015(Nijhout et al. , 2018;;Best et al. 2009;Nijhout and Reed 2014;Nijhout et al. 2017Nijhout et al. , 2019) ) among others, have shown that homeostasis is an important phenomenon in biochemical reaction networks.In a biochemical network, each node represents the concentration of a chemical substrate and each arrow denotes a chemical interaction between the molecules at the head and tail of the arrow.In an input-output network formulation, one node is designated as the input node ι and another is designated as the output node o.The modeling assumes that some external stimuli (represented by an input parameter, or simply an input, I) affects the network dynamics only at the input node, and the end result of computation by the network dynamics is the value of the output node.In this setting, there is a canonical choice for the smooth observable φ: the coordinate function of the output node.
Motivated by these examples, Wang et al. (2021) introduced the notion of 'abstract input-output network' and devised a scheme for the classification of 'homeostasis types' in such networks.The notion of homeostasis type of a network makes precise the idea that homeostasis may be caused by different 'mechanisms' in that network.The results of Wang et al. (2021) apply to the case of single-input single-output networks where the external stimuli can only affect one input node via a single scalar input.
Even though single-input single-output (SISO) networks are quite popular in many engineering domains (Ma et al. 2009;Ang and McMillen 2013;Araujo and Liotta 2018;Bhattacharya et al. 2022), the single input node and single input assumptions seem unrealistic in biology, as disturbances that arise are typically very complex and do not have a single well-defined entry point (Gupta and Khammash 2022).As far as the input is concerned, there are two possible ways to extend the work of Wang et al. (2021) in order to include more complex situations: (1) Multiple input nodes.A single input affects more than one of several input nodes.
(2) Multiple inputs.Several inputs affect more than one of several input nodes.
Regarding (1), Madeira and Antoneli (2022) extended the classification of Wang et al. (2021) to the setting of multiple input nodes and used this extended theory to completely work out the homeostasis types of a representative model for bacterial chemotaxis (Clausznitzer et al. 2010;Tindall et al. 2008).As for (2), an interesting biologically relevant example is the regulation of extracellular dopamine (eDA) in response to variation in the activities of the enzyme tyrosine hydroxylase (TH) and the dopamine transporters (DAT) (Best et al. 2009;Golubitsky and Stewart 2018).Another biologically relevant example of the second situation is the mathematical model of Granjon et al. (2017) for the physiological co-regulation of calcium (Ca) and phosphate (PO 4 ) in the rat.We will discuss this model in details in Sect.3.
Let us recall the results of Madeira and Antoneli (2022).The main discovery of Wang et al. (2021) is that, in a given abstract input-output network, there is a finite number of 'distinct mechanisms' that may cause homeostasis, i.e. may force the derivative of the input-output function to vanish (at a fixed value of the input).These 'distinct mechanisms', called 'homeostasis types', bijectively correspond to a specific subnetworks of the abstract input-output network.These subnetworks, called homeostasis subnetworks, can be characterized in purely topological terms.In the single-input single-output theory of Wang et al. (2021), the homeostasis subnetworks can be divided into two classes: structural and appendage.The structural subnetworks correspond to feedforward mechanisms and the appendage subnetworks correspond to feedback mechanisms.These are called, respectively, 'opposer modules' and 'balancer modules' in the control theoretic literature (Araujo and Liotta 2018).In Madeira and Antoneli (2022) the single-input single-output theory is generalized, for the first time, to account for multiple input nodes.The main result is that everything from the singleinput case generalizes, except that there is a new class of homeostasis subneworks, called input counterweight: every abstract multiple input node network has exactly one input counterweight subnetwork, which can be topologically characterized, as well.
In this paper we further build on the theory of Wang et al. (2021); Madeira and Antoneli (2022) to completely solve the problem of classifying the homeostasis types for input-output networks with multiple input nodes and multiple inputs (Sect.2.3).More precisely, given an input-output network with multiple inputs we show that one can consider each parameter at a time.Thus, effectively reducing the problem of classification of homeostasis types to the single input case (still with multiple input nodes), that has been completely solved in Madeira and Antoneli (2022).Afterwards, we show how to combine these partial classifications for the single input cases into an algorithm that provides the full classification on the multiple inputs setting (Sect.2.8).
The first issue to consider in the multiple inputs case is to obtain an analogue of the 'determinant formula' for the gradient of the inputs-output function (cf.Golubitsky and Stewart 2017b;Golubitsky et al. 2020;Wang et al. 2021;Madeira and Antoneli 2022).Versions of this 'determinant formula' have been obtained by several authors: Ma et al. (2009), Golubitsky and Wang (2020) for three-node networks and Araujo and Liotta (2018), Aoki et al. (2019) for arbitrary networks under the name 'RPA equation'.In Sect.2.3 we introduce a definition of a multiple-input single-output network and prove a multivariate generalization of 'determinant formula' (Lemma 2.1) mentioned above.A similar result to our 'determinant formula' has been obtained in Tang and McMillen (2016).
The main result of this paper allows us to completely classify homeostasis subnetworks of multiple inputs 'core' network (Sects.2.4, 4.1).In the multiple inputs setting two new features arise.The first is related to being able to consider one parameter at a time.More specifically, for each scalar input I M we define a unique subnetwork, called I M -specialized subnetwork, that contains all the homeostasis subnetworks associated to the input I M (Sect.2.5).When we have all the homeostasis subneworks of the specialized subnetworks we can proceed by considering two cases: (i) a pleiotropic subnetwork if it appears in all I M -specialized subnetworks, (ii) a coincidental subnetwork otherwise.Then we show that a pleiotropic subnetwork can only be of structural or appendage classes (Sects.2.6, 4.2), whereas a coincidental subnetwork can be of structural, appendage or input counterweight class (Sect.2.7).Furthermore, all our results hold generically in the setting of influence networks, which includes all ODE models used in biology (see Sect. 2.2).
The second new feature of multiparameter setting is related to the occurrence of overlapping between coincidental subnetworks contained in distinct I M -specialized subnetworks.These non-trivial interactions between homeostasis subnetworks in multiple inputs networks leads to the appearence of homeostasis mode interaction or higher codimension homeostasis.The notion of mode interaction is familiar in bifurcation theory.In a steady-state bifurcation the eigenvectors of the linearized equation corresponding to simple eigenvalues are called modes.A mode whose eigenvalues lie on the imaginary axis is said to be critical.Generically, it is expected that an one-parameter system has only one critical mode.However, in systems with more than one parameter, one expects multiple critical modes.The steady-state bifurcations that may arise in nonlinear systems near a (multi-) parameter value at which there are multiple critical modes are thought of as resulting from a nonlinear interaction of several critical modes.This process is called mode interaction and the (multi-) parameter values at which there are multiple critical modes are called higher codimension bifurcation points.We assert that there is an analogue process in the context of infinitesimal homeostasis (Ex.2.16).Duncan et al. (2023) investigates the appearance of codimension-two homeostasis mode interaction in the different setting of infinitesimal homeostasis with a single input parameter.
Structure of the Paper.In Sect. 2 we give the definitions and state the main results of the paper.We discuss some simple abstract examples to illustrate the definition and results.In Sect. 3 we apply our results to a 'real biological system', a mathematical model of calcium and phosphate metabolism proposed in Granjon et al. (2017).In Sect. 4 we give the proofs of all our results.Finally, in Sect. 5 we briefly discuss our results in the context of the theory of infinitesimal homeostasis and conclude the paper with an outlook for future research.

Homeostasis in multiple inputs networks
In this section we state the main results of the paper and provide the necessary definitions.

Dynamical theory of infinitesimal homeostasis
Golubitsky and Stewart proposed a mathematical method for the study of homeostasis based on dynamical systems theory (Golubitsky andStewart 2017b, 2018) [see the review (Golubitsky et al. 2020)].In this framework, one consider a system of differential equations where X = (x 1 , . . ., x n ) ∈ R n is the vector of state variables and I = (I 1 , . . ., I N ) ∈ R N is the vector of input parameters.Suppose that (X * , I * ) is a linearly stable equilibrium of (2.1).By the implicit function theorem, there is a function X (I) defined in a neighborhood of I * such that X (I * ) = X * and F( X (I), I) ≡ 0. See Jahedi et al. (2022) for results on the generic existence and robustness of X (I).
A smooth function φ : R n → R is called an observable.Define the input-output function z : R N → R associated to φ and X as z(I) = φ( X (I)).The input-output function allows one to formulate several definitions that capture the notion of homeostasis (see Ma et al. 2009;Ang and McMillen 2013;Tang and McMillen 2016;Golubitsky andStewart 2017b, 2018).
Definition 2.1 Let z(I) be the input-output function associated to a system of differential Eqs.(2.1).We say that z(I) exhibits (a) Perfect Homeostasis on an open set ⊆ dom(z) if That is, z stays within the range z(I s ) ± δ over .(c) Infinitesimal Homeostasis at the point It is clear that perfect homeostasis implies near-perfect homeostasis, but the converse does not hold.Inspired by Nijhout, Reed, Nijhout et al. (2004); Best et al. (2009); Golubitsky andStewart (2017b, 2018) introduced the notion of infinitesimal homeostasis that is intermediate between perfect and near-perfect homeostasis.It is obvious that perfect homeostasis implies infinitesimal homeostasis.On the other hand, it follows from Taylor's theorem that infinitesimal homeostasis implies near-perfect homeostasis in a neighborhood of I c , see Golubitsky and Stewart (2022).It is easy to see that the converse to both implications is not generally valid (see Reed et al. 2017).
The notion of infinitesimal homeostasis allows one to apply the tools from singularity theory.For instance, by considering higher degeneracy conditions, in addition to (2.4), one is lead to distinct forms of infinitesimal homeostasis that can be classified by elementary catastrophe theory (see Golubitsky andStewart 2017b, 2018 for details).Finally, when combined with coupled systems theory (Golubitsky and Stewart 2006) the formalism of Golubitsky and Stewart (2017b), Golubitsky and Stewart (2018), Golubitsky et al. (2020) becomes very effective in the analysis of model equations.

Networks and dynamical systems
Before defining the appropriate class of dynamical systems (in Sect.2.3) for our results we will briefly discuss the relation between networks and dynamics.
A large portion of the literature on network dynamical systems modeling seems to suggest that there is a unique way to associate a system of differential equations to a given directed graph G.With some rare exceptions, e.g.Bick et al. (2023), a precise definition of what is a "network dynamical system" is completely overlooked and the pairwise interaction interpretation (see below) is assumed without further justification.In fact, there are at least, two possible ways to attach a (class of) dynamical system(s) to a directed graph G.In order to discuss the distinction between these two possibilities we need precise definitions.Let G be a directed graph with k nodes.
(1) Pairwise interaction interpretation.In this interpretation the dynamics is encoded by a weighted adjacency matrix A compatible with the directed graph G.That is, A ji = 0 if and only if there is a link from node j to node i.For simplicity, suppose that each node of G represent an identical dynamical system of the form ẋi = F(x i ), where x i ∈ R n is the state vector of node i and F : R n → R n is a smooth function that describes the internal dynamics of node i.The interaction between nodes is given by a smooth pairwise coupling function G : R n × R n → R n , also called point-to-point coupling (Golubitsky and Stewart 2002).The tuple (A, F, G) defines a network dynamical system through the set of differential equations The dynamics is determined by the evolution of the joint state of all nodes (x 1 , . . ., x k ) through (2.5).It is important to note that network dynamical systems described via (2.5) have only additive interactions.Specifically, the interactions are in general nonlinear in the state variables x i , but linear in the coupling weights A ji Aguiar and Dias (2018).By letting F and G vary over all smooth functions on the state space R n , for a given family of weighted adjacency matrix A defining the same directed graph G, one obtains a space of vector fields X ∞ A (R n ).While setup (2.5) is arguably one of the most commonly used formulations of network-based modeling, it imposes a severe restriction that might not hold in real-world systems (Bick et al. 2023) (Golubitsky andStewart 2006, 2022).Now, despite node dependencies being captured by a graph, this does not exclude the possibility of nonlinear interactions involving three or more nodes, called higher-order interactions.By letting F and H vary over all smooth functions on the state space R n one obtains a space of vector fields X ∞ G (R n ), this is the space of all dynamical systems (2.6) that are compatible with the network structure G.The main goal of the groupoid formalism is to study the space X ∞ G (R n ), rather than considering (2.6) for a specific F and H .This yields insights on how the 'generic' dynamical behavior of such a system depends on the imposed network structure encoded by G.
It is easy to see that the space of vector fields X ∞ G (R n ) contains the space of vector fields X ∞ A (R n ), for all families of weighted adjacency matrices A defining the directed graph G.In fact, X ∞ G (R n ) is the largest space of vector fields that can be attached to a directed graph G.Likewise, X ∞ A (R n ), for a fixed A, is the smallest space of vector fields that can be attached to a family of weighted adjacency matrix A. Remarkably, it seems that in order to generalize (2.5) to include higher order interactions only up to a fixed level, directed graphs are not enough to capture all the interaction relations.It is necessary to consider 'higher dimensional' generalizations of graphs, such as hypergraphs or simplicial complexes (Bick et al. 2023).
In this paper we adopt the influence network interpretation.As mentioned before, this approach emphasizes the 'generic' properties and how they are affected by the network structure encoded by G.This leads to the notion of model-independent approach, which is thoroughly explored in the book Golubitsky and Stewart (2022).The modelindependent approach contrasts with the more common model-dependent one, where specific model equations are solved, usually numerically.Both approaches have advantages and disadvantages and should not be seen as competing against each other, on the contrary, they complement each other.Specific models are useful in connection with experimental tests, when the equations are obtained from precise laws.In a modelindependent approach the exact equations may not be known, nevertheless one still can predict, from known properties of the equations, what types of behaviors are expected and which ones are forbidden.

Multiple-input single-output networks
A multiple-input single-output (MISO) network, or simply a multiple inputs network, is a network G with n distinguished input nodes ι = {ι 1 , ι 2 , . . ., ι n }, all of them associated to at least one input parameter I M , M = 1, . . ., N , one distinguished output node o, and r regulatory nodes ρ = {ρ 1 , . . ., ρ r }.The associated system of differential equations have the form is the vector of state variables associated to the network nodes.
We write a vector field associated with the system (2.7) as and call it an admissible vector field for the network G. Let f j,x i denote the partial derivative of the j th node function f j with respect to the i th node variable x i .We make the following assumptions about the vector field F throughout: (a) The vector field F is smooth and has a linearly stable equilibrium at (X * , I * ).
Therefore, by the implicit function theorem, there is a function The partial derivative f j,x i can be non-zero only if the network G has an arrow i → j, otherwise f j,x i ≡ 0. (c) Only the input node coordinate functions f ι k depend on at least one of the components of the vector of input parameters I and the partial derivative of f ι k ,I M generically satisfies for some M = 1, . . ., N .

Remark 2.2
In this paper we explicitly exclude the possibility that the output node is one of the input nodes.This assumption is included purely for the sake of convenience.
In fact, all the results should be valid in this case, but then all the theorems and proofs should be properly adapted to take this particular case into account.This possibility, in the case of single input networks will be considered with great detail in another publication (Antoneli et al. 2023).♦ The network structure provides a distinguished class of observables for an admissible system, namely, the state variables.In the particular case of an input-output network the observable of interest is given by the output state variable x o .
Definition 2.3 Let G be a multiple inputs network and F be a family of admissible vector fields with an equilibrium point Infinitesimal homeostasis in a multiple inputs network is given by the critical points of x o (I), namely, the zeros of the gradient vector By a straightforward application of Cramer's rule Wang et al. (2021) obtained a determinant formula for the derivative of the input-output function in the single-input single-output case.Madeira and Antoneli (2022) generalized the determinant formula of Wang et al. (2021) to the case of multiple input nodes networks.In the following we further generalize the determinant formula of Madeira and Antoneli (2022) to the case of multiple inputs networks.
Let J be the (n + r + 1) × (n + r + 1) Jacobian matrix of an admissible vector field Here all partial derivatives f ,x j are evaluated at X (I), I .

Lemma 2.1 Let x o (I) be input-output function of a multiple inputs network. The partial derivative of x o (I) with respect to the M-th parameter
where det(J ) and det H M are evaluated at the equilibrium point X (I).Hence, Moreover, I 0 is a point of infinitesimal homeostasis if and only if as a function of I evaluated at I 0 .
Proof Implicit differentiation of the equation F( X (I), I) = 0 with respect to I yields the linear system Since X (I) is assumed to be a linearly stable equilibrium, it follows that det(J ) = 0. On applying Cramer's rule to (2.15) we can solve for ∂ x o ∂I M (I) obtaining (2.12).Applying (2.12) to (2.9), we obtain Eq. (2.13).
Remark 2.4 An explicit expression for det H M can be obtained by expanding it with respect to the last column and the ι m -th row: where H ι m is obtained from H by removing the last column and the ι m -th row.When there is a single input, i.e.N = 1, the gradient ∇x o reduces to ordinary derivative x o and (2.13) gives the formula for x o obtained in Madeira and Antoneli (2022).When there is a single input and a single input node, N = n = 1, there is only one matrix H ι m = H , called the homeostasis matrix and (2.13) gives the corresponding formula for x o obtained in Wang et al. (2021).♦

Core networks
Wang et al. ( 2021) introduced a fundamental construction for the study of homeostasis in input-output networks, called 'core subnetwork'.Madeira and Antoneli (2022) extended this construction to the case of input-output networks with multiple input nodes and single input parameter.Here we extend the arguments of Wang et al. (2021), Madeira and Antoneli (2022) to the case of multiple inputs networks.Let (X * , I * ) be a linearly stable equilibrium of the G-admissible ODE system (2.1).Then (X * , I * ) it satisfies the system of equations Partition the nodes of the network G as follows: (i) input nodes (whose dynamics explicitly depends on at least one input), (ii) the output node and (iii) the regulatory nodes, that can be classified into three types depending if they are upstream from the output node or/and downstream from at least one input node.More precisely, the set of regulatory nodes may be partitioned as: (a) Those nodes σ that are both upstream from o and downstream from at least one input node ι m , (b) Those nodes d that are not downstream from any input node ι m , (c) Those nodes u which are downstream from at least one input node ι m , but not upstream from o. Proof Follows from Theorem 4.2.
Theorem 2.2 allows one to assume, without loss of generality, that G is a core network, that is G = G c , as far as infinitesimal homeostasis is concerned.

Structure of infinitesimal homoestasis
In this subsection, unless explicitly stated, we assume that G is a core multiple inputs network with input nodes ι 1 , . . ., ι n and inputs I 1 , . . ., I N .
By Lemma 2.1 a network G exhibits infinitesimal homeostasis at a point I 0 whenever the vector-valued function (when evaluated at ( X (I), I)) vanishes at I 0 : (2.18) here det H M are the determinants of the I M -generalized homeostasis matrices.
In order to analyze and simplify these determinants let us introduce some terminology.A multivariate vector-valued polynomial, or, simply a polynomial mapping is a mapping P : R k → R k with polynomial components.That is, if P is a polynomial mapping, there exist multivariate polynomials P 1 , P 2 , . . ., P k : R k → R such that we can write P as We say that P is irreducible if and only if each component P j is irreducible.Suppose there is a multivariate polynomial function p : R k → R that is common factor to all components P j .Then we can factor p from the polynomial vector P as We say that p is a scalar factor of P.
Recall that the nonzero entries of the I M -generalized homeostasis matrices H M are the partial derivatives f j,x i and f j,I M .In particular, det H M is a homogeneous polynomial function of degree (n + r + 1) in the partial derivatives f j,x i and f j,I M .Hence, the vector-valued function h is a (formal) polynomial mapping on the 'variables' f j,x i and f j,I M .The scalar-valued function h(I) (depending on the vector of input parameters I) is obtained by evaluating the partial derivatives f j,x i and f j,I M at X (I).
Let us motivate the next definition with a simple observation.In the multiparameter setting, even a core network may have nodes that are not affected by all inputs.For example, consider the 5-node multiple inputs network G shown in Fig. 2. In this figure inputs are highlighted by dotted circles of distinct colors.The arrows from each input to distinct input nodes are of the same color (the same color of the corresponding dotted circle).The network shown in (a) has three input nodes, ι 1 , ι 2 and ι 3 , and two Fig. 2 A 2-parameter 5-four node network (a) and its I M -specialized subnetworks (b) and (c).Here, inputs are highlighted by dotted circles of distinct colors.The arrows from each input to distinct input nodes are of the same color (the same color of the corresponding dotted circle).The 2-parameter network (a), has 5 nodes: three input nodes ι 1 , ι 2 , ι 3 , one output node o and one regulatory node σ .The two I M -specialized subnetworks of (a) are always single-input multiple input node networks (see Definition 2.6).The I 1specialized subnetwork (b) has 1 input node (ι 1 ).Notice that ι 2 becomes a regulatory node for this network.The I 2 -specialized subnetwork (c) has 2 input nodes (ι 2 , ι 3 ) inputs, I 1 (blue) and I 2 (red).Although, G is a core network, the input node ι 1 is not affected by the parameter I 2 , and the node ι 3 is not affected by parameter I 1 .To overcome this difficulty, we define the 'specialized networks' relative to a single input parameter.
Definition 2.6 Let G be a core multiple inputs network with inputs I 1 , . . ., I N .The I M -specialized subnetwork G I M is defined as the (single input) input-output subnetwork of G consisting of all the input nodes that receive the input I M , all the regulatory nodes that are downstream from those input nodes and the output node.The arrows of G I M are the arrows of G between the nodes of G I M .The subnetwork ♦ The specialized subnetwork G I M can be considered as a multiple input node network with single input I M , as studied in Madeira and Antoneli (2022), by 'forgetting' the effect of the other parameters and reducing to the core network using the core reduction theorem of (Madeira and Antoneli 2022, Thm. 3.2).The input nodes of G I M are exactly the input nodes of G that are affected by the parameter I M .Definition 2.7 Let G be a core multiple inputs network with inputs I 1 , . . ., I N .The I M -homeostasis matrix H I M associated to G I M , is the generalized homeostasis matrix of the multiple input-output network G I M (see Madeira and Antoneli (2022, Sec 2.3)).Similarly, the corresponding I M -vestigial subnetwork D I M has an associated jacobian matrix J D I M (see Madeira and Antoneli (2022, Sect. 3.2)).To simplify the notation, in the case where the I M -vestigial subnetwork is empty (D I M = ∅), we write det(J D I M ) ≡ 1. ♦ Now we need to specify the set ⊆ R N of allowable parameter values.This set depends on the admissible vector field F and the type of model being considered.For instance, in a biochemical network the set is the positive orthant of some R N .The subset of non-singular parameters of F on is defined as (2.19) where J = (D F) ( X (I),I) is the jacobian.The set J also depends on the vector field F and, generically, is an open dense subset of .

Lemma 2.3
For each M = 1, . . ., N , we have: Lemma 2.3 allows us to further simplify the components of h, by considering the determinants det H I M .This reduce to the situation already studied in Madeira and Antoneli (2022).
Definition 2.8 The vector determinant associated to an input-output network is the vector-valued function defined by where det H I M , M = 1, . . ., N , is the determinant of generalized homeostasis matrix of the I M -specialized subnetwork G I M .The vector-valued function h can be considered as a (formal) polynomial mapping on the 'variables' f j,x i and f j,I M .♦ Proposition 2.4 The vector-valued functions ∇x o , h and h, defined on → R N , have the same set of zeros on J .
Proof The first equality follows from Lemma 2.1 and the second equality follows from Lemma 2.3.
The König-Frobenius theorem (Schneider 1977;Brualdi and Cvetkoić 2009) (see also Wang et al. (2021); Madeira and Antoneli (2022)) imply that the components of the polynomial mapping h can be factorized as the product of the determinants of the irreducible diagonal blocks of each H I M (defined up to row and column permutations).An irreducible block B of some H I M is called a homeostasis block.We can further collect the factors that are common irreducible diagonal blocks of all matrices H I M and bring them to the front as scalar factors.Then we can write Therefore, we can split the problem of classifying homeostasis types supported by G into two cases according to whether the components of h have a common scalar factor or not.
Definition 2.9 Let G be a core multiple inputs network and consider its vector determinant h as in (2.22).A homeostasis block corresponding to scalar factor det(B i ) The other homeostasis blocks of G are called coincidental.♦ Remark 2.10 In genetics, pleiotropy refers to the phenomenon when a single locus affects multiple traits (Stearns 2010).Here, we employed the term pleiotropic homeostasis referring to the fact that the nullification of one single homeostasis block leads to the annulment of the whole homeostasis vector h.♦ Recall that the homeostasis types of a single parameter input-output network G are given in terms of the factors of h = det(H ) (see Wang et al. (2021); Madeira and Antoneli (2022)) where each irreducible block B j can be of three types, called appendage, structural and counterweight.There is only one counterweight block defined as the only irreducible block that contains all the partial derivatives of f with respect to the input I.
Infinitesimal homeostasis of type B j occurs if det(B j ) = 0 and det(B i ) = 0 for all i = j.This is generic when there is only one input parameter.Now, suppose that G has N inputs affecting n input nodes.Generically, it is expected that N irreducible factors of h can simultaneously vanish at a fixed input value.

Definition 2.11
Let G be a multiple inputs core network.Let B be a homeostasis block of size .We say that the homeostasis class of B is

Pleiotropic homeostasis types
In this section we classify the pleiotropic sub-blocks of a multiple inputs core network.(ii) This follows immediately from the fact that the classification is based on the number of self-coupling entries.
In Remark 2.4 we observed that when the network G has a single parameter the theory developed here reduces to the situation considered in Madeira and Antoneli (2022).The next result shows that when the multiple inuts network G has a single input node then the theory essentially reduces to the case where there is only one input parameter, considered in Wang et al. (2021).This is an extreme case where only pleiotropic homeostasis occurs.
Proposition 2.6 Suppose the core multiple inputs network G has only one input node ι and multiple inputs I 1 , . . ., I N .Then, we have: Therefore, Proposition 2.6 implies that the classification of homeostasis types of a multiple inputs network with a single input node is exactly the same as that of a single input parameter single input node network.However, it is not true that we expect to see the same 'homeostasis phenomena' in both networks.In order to see this, suppose that G has a single input node, but multiple inputs (I 1 , . . ., I N ) affecting the input node.Consequently, the function h, which is the same polynomial on the partial derivatives f j,x i , becomes a multivariate function when evaluated at X (I 1 , . . ., I N ).Generically, it is expected that N irreducible factors of (2.23) can simultaneously vanish at a fixed value (I 0 1 , . . ., I 0 N ).Returning to the general case, the next result gives a complete topological classification of the homeostasis subnetworks corresponding to the pleiotropic homeostasis types.
Theorem 2.7 Let G be a multiple inputs core network, B be a pleiotropic block of G and K B be the corresponding homeostasis subnetwork to B in G.
(i) Then K B is either an appendage or structural subnetwork of all I M -specialized subnetworks G I M .(ii) More precisely, K B is an appendage (respect.structural) subnetwork if and only if it is a H I M -appendage (respect.H I M -structural) subnetwork, for some (and hence for all) M = 1, . . ., N .
Proof The result follows from Theorems 4.3 and 4.4 for the appendage case and Theorems 4.8 and 4.9 for the structural case.See Sect.4.2 for precise characterization of each type of subnetwork.

Coincidental homeostasis types
The occurrence of coincidental homeostasis reflects the fact that the mechanism leading to homeostasis is not the same with respect to all the inputs.The coincidental homeostasis types are given by all the possible combinations of coincidental blocks.
A coincidental type can be of structural class, appendage class or input counterweight class.Thus, depending on the network, a coincidental type can have the form of a mix of these three classes.
In the simplest case the determinant of a coincidental block can be an entry of some H I M of the form f x ιm ,I M .Since, by assumption, these entries cannot vanish it may happen that some coincidental types do not yield infinitesimal homeostasis.On one hand, as shown in Proposition 2.6, networks that have only one input node, only support pleiotropic homeostasis.On the other hand, it is easy to find networks that support only coincidental homeostasis (see Sect. 2.9).
Still, one may wonder whether there is a multiple inputs core network that do not support either pleiotroic nor coincidental homeostasis, that is, it does not support infinitesimal homeostasis.The next proposition shows that this cannot happen, namely, any multiple inputs core network always support at least one type of homeostasis.

Proposition 2.8 A multiple inputs core network G always supports infinitesimal homeostasis.
Proof In order to prove the proposition, consider all input nodes ι m , such there is an ι m -simple node σ (σ = ι m ) that receives an arrow from ι m and such that f σ,x ιm = 0 (Haldane homeostasis).Let G m be the core subnetwork between the input node ι m and output node o (Madeira and Antoneli 2022, Def. 2.13).By definition, G m is a single input-output network.Then, by the characterization of structural homeostasis in networks with single one input node (see Wang et al. 2021), the homeostasis determinant of G m vanish, i.e., if H c ι m is the homeostasis matrix of G m , then det H c ι m = 0.This fact together with (Madeira and Antoneli 2022, Eq. 3.39) implies that det H I M = 0, for all M = 1, . . ., N .To conclude the argument, we claim that this construction does not force det(J ) = 0, generically.This follows from that fact that one of the terms that appear in the expression of the Jacobian determinant is the product of all the self-couplings of nodes of G.As the construction above does not assume that the self-couplings to be equal to zero, the claim holds.
Finally, we state below a sufficient condition for a core multiple inputs network to support coincidental homeostasis.The network of Fig. 4b shows that the condition given below is sufficient, but not necessary, for a core multiple inputs network to support coincidental homeostasis.

Proposition 2.9 Let G be a multiple inputs core network. If none of the input nodes of G is an absolutely super-simple node (see Definition 4.1), then G supports at least one coincidental homeostasis type.
Proof Note that an input node is an absolutely super-simple node if and only if it is the first absolutely super-simple node of the network.Hence, the hypothesis that none of the input nodes of G is an absolutely super-simple node is equivalent to: the first absolutely super-simple node of G is not an input node.First, suppose that pleiotropic homeostasis does not occurs in G. Then by Proposition 2.8 it follows that coincidental homeostasis must occur in G. Now suppose that pleiotropic homeostasis does occur in G. Clearly, as the input nodes are not absolutely appendage, by Theorem 2.7, they do not belong to a pleiotropic-appendage subnetwork.In addition, none of the input nodes belong to any absolutely super-simple structural subnetwork, as the input nodes are not absolutely super-simple and cannot be between two absolutely super-simple Remark 2.14 As previously explained, Proposition 2.9 does not provide us with a necessary condition for the occurrence of coincidental homeostasis.It is then a open problem to determine necessary and sufficient conditions for the occurrence of coincidental homeostasis.As it will be exemplified in Sect.2.9, the interaction between the different subnetworks makes this a non-trivial problem.♦

Algorithm to determine all homeostasis types
Using the results obtained here, together with Wang et al. (2021); Madeira and Antoneli (2022), one can write down a general algorithm to find all homeostasis types of a given multiple inputs core network G.
Step 1 For each parameter I M with M = 1, . . ., N determine the I M -specialized subnetwork G I M as in Definition 2.6.
Step 2 Since each G I M is a single parameter multiple input-output network one can apply the algorithm of (Madeira and Antoneli 2022, Sect. 2.6) to determine all the homeostasis subnetworks K of each G I M .
Step 3 Determine the homeostasis subnetworks that are common to all G I M .The output of this step is the list of pleiotropic homeostasis types.
Step 4 Determine the N -tuples formed by combinations of the remaining homeostasis subnetworks.The output of this step is the list of coincidental homeostasis types.

Some simple examples
In this section we present some examples of small networks exhibiting an astonishing array of phenomena that can arise in the multiparameter setting.
For instance, there are networks that do support both pleiotropic and coincidental homeostasis (e.g.Example 2.16) and networks that support only one of each type: Example 2.15 supports only pleiotropic homeostasis.Whereas Example (...) supports only coincidental homeostasis.
Moreover, we will see that when there is a proper coincidental type (i.e.all determinants can vanish) non-trivial interaction between the homeostasis subnetworks can occur.
Let us start with the minimal non-trivial examples.There are some constraints to get non-trivial networks: (a) There must be at least two inputs.(b) There must be at least two input nodes (see Proposition 2.6).(c) The output node is distinct from the input nodes (see Remark 2.2).Therefore, the network must have at least three nodes, two input nodes ι 1 , ι 2 and the output node o.In addition, let us make the simplifying assumption that each input parameter affects exactly one input node.Granted these conditions, it is not difficult to show that, up to core-equivalence and input-relabeling, there are exactly 8 networks, shown in Fig. 3. Two core networks are core-equivalent if they have identical vector determinants, up to permutation (cf.Wang et al. 2021, Def. 1.9).
Example 2.15 Consider the eight 3-node core networks shown in Fig. 3.We begin by determining the specialized subnetworks.Since the networks have one input per input node the specialized subnetwork is a single parameter single input 3-node network.Then we can use the results from Golubitsky and Wang (2020) to write the vector determinant.We also note that, by definition, f ι 1 ,I 1 and f ι 2 ,I 2 cannot vanish.Network (a).In this example we have Therefore, the vector determinant is given by Hence, the network only supports coincidental homeostasis, given by the simultaneous vanishing of the irreducible (structural) factors f o,x ι 1 and f o,x ι 2 .

Network (b).
In this example we have In the second specialized network we omitted the arrow o → ι 2 , since it does not affect the vector determinant.Therefore, the vector determinant is given by Hence, the network does not support pleiotropic homeostasis.There are two possibilities for the occurrence of coincidental homeostasis (1) f o,x ι 1 (structural) and f o,x ι 2 (structural), (2) f ι 2 ,x ι 2 (appendage) and f o,x ι 2 (structural).Now, the jacobian determinant of the network is Hence, the occurrence of homeostasis in case (2) necessarily forces det(J ) = 0 (see Remark 2.17).Therefore, coincidental homeostasis only can occur in case (1).

Network (c).
In this example we have Therefore, the vector determinant is given by Now, the jacobian determinant of the network is Hence, the network only supports coincidental homeostasis.Moreover, from the 4 possible combinations of factors of h the only that does not force det(J ) = 0 is given by f o,x ι 1 and f o,x ι 2 (see Remark 2.17).

Network (d).
In this example we have In both specialized networks we omitted the arrow o → ι 2 since it does not affect the vector determinant.Therefore, the vector determinant is given by Hence, the network only supports coincidental homeostasis, given by the simultaneous vanishing of the two irreducible (structural) factors.
Network (e).In this example we have In the second specialized network we omitted the arrow o → ι 2 , since it does not affect the vector determinant.Therefore, the vector determinant is given by Hence, the network does support pleiotropic homeostasis.In order for coincidental homeostasis to occur, both irreducible factors below must vanish simultaneously Now the jacobian determinant of the network is Thus, the vanishing of the two factors force det(J ) to vanish and so coincidental homeostasis can not occur in this network (see Remark 2.17).

Network (f).
In this example we have Therefore, the vector determinant is given by Hence, the network does not support pleiotropic homeostasis.In order for coincidental homeostasis to occur, both irreducible factors below must vanish simultaneously Note that both factors are of structural homeostasis type.The vanishing of f o,x ι 2 reduces the first factor to f o,x ι 1 f ι 2 ,x ι 2 .The jacobian determinant of the network is Hence, the condition f ι 2 ,x ι 2 = 0 forces det(J ) = 0, whereas f o,x ι 1 = 0 does not (see Remark 2.17).Therefore, the vanishing of the factor f o,x ι 1 0 is the only possibility for the occurrence of coincidental homeostasis.

Network (g).
In this example we have Therefore, the vector determinant is given by Hence, the network does not support coincidental homeostasis.In order for pleiotropic homeostasis to occur the irreducible (structural) factor f o,x ι 2 must vanish.Here, the obstruction to the occurrence of coincidental homeostasis is due to the fact that the second component of ĥ consists only of the non-vanishing factor f ι 2 ,I 2 .Network (h).In this example we have In the first specialized network we omitted the arrow ι 2 → ι 1 , since it does not affect the vector determinant.Therefore, the vector determinant is given by Hence, pleiotropic homeostasis occurs when the irreducible (structural) factor f o,x ι 2 vanishes.Coincidental homeostasis is given by the simultaneous vanishing of the irreducible factors f ι 2 ,x ι 1 and f ι 1 ,x ι 1 .♦ Next, we consider an example obtained from the network shown in Fig. 3g by adding the influence from input I 2 on node ι 1 .The main difference from the networks in Example 2.15 is the appearance of an input counterweight factor.This illustrates the effect of having more than one input node affected by the same input (see Fig. 4).Note that this is core equivalent to adding the influence from input I 2 on node ι 1 on the network shown in Fig. 3h.

Example 2.16
Consider the 2-parameter 3-node core network shown in Fig. 4. The specialized subnetworks are: Therefore, the vector determinant is given by Pleiotropic homeostasis can occur by the vanishing of the irreducible (structural) factor f o,x ι 2 .In order for coincidental homeostasis to occur both irreducible factors (the first is structural and the second is input counterweight) below must vanish simultaneously Occurrence of coincidental homeostasis requires the simultaneous vanishing of f ι 2 ,x ι 1 and f ι 1 ,x ι 1 .Now, the jacobian determinant of the network is Hence the vanishing of the two factors above forces det(J ) = 0. Thus, coincidental homeostasis cannot occur in this network (see Remark 2.17).♦ Remark 2.17 In three of the eight networks of Example 2.15, networks (b), (e), (f), and the network of Example 2.16, we faced a situation where the vanishing of certain factors, that could cause coincidental homeostasis, forced the vanishing of the jacobian determinant of the network.In other words, the homeostasis point occurs at the same value of the input parameters that makes the family of equilibria to undergo a steady-state bifurcation.Strictly speaking, this means that the situation mentioned above cannot be considered as a 'proper' infinitesimal homeostasis and thus we have excluded these cases.Indeed, the definition of infinitesimal homeostasis (Definition 2.1) excludes the simultaneous occurrence of homeostasis and steadystate bifurcation.However, if one considers extending the definition of homeostasis to include such cases (see Duncan et al. (2018); Duncan and Golubitsky (2019) for some advances in this direction) then one may get a much richer variety of phenomena.♦

Application of the theory to a model of calcium and phosphate homeostasis
In this section, we shall apply our theory to a 'real biological system', a mathematical model for the metabolic regulation of calcium and phosphate proposed by Granjon et al. (2017).Our aim here is to exemplify an application of the results of this paper to mathematical model from the literature, rather than exploit the details of calcium metabolism under physiological and pathological conditions.
Calcium is an essential metal ion that takes part in many signaling pathways and biochemical processes, including bone metabolism.Hence, its extracellular concentration must be tightly regulated (Blaine et al. 2015).Importantly, the regulation of extracellular calcium concentration is coupled to phosphate homeostasis (Blaine et al. 2015), which is an anion essential to human body.An explanation of the many pathways involved in calcium and phosphate metabolism is beyond the scope of this paper, and the interested reader is referred to Blaine et al. (2015) and Melmed et al. (2015, Ch 29).
The model used here is adapted from Granjon et al. ( 2017) and can be described by a network with 7 nodes (two input nodes, one output nodes and four regulatory nodes) and two inputs (see Fig. 5).The two inputs represent (1) the calcium intake (I 1 ) and ( 2) the phosphate intake (I 2 ).This is a quite natural choice, since one can assume that the intake of these ions depend on the availability of food, which can be reasonably variable.These inputs directly affect, respectively, the dynamics of the intestinal concentration of free calcium (x ι 1 ) and of free phosphate (x ι 2 ).Hence, these concentrations correspond to the two input nodes.The output node corresponds to the extracellular concentration of calcium x o .Since calcium and phosphate metabolism is regulated by a complex network of hormones and signaling molecules, we include here the components that were highlighted by Granjon et al. (2017).The four regulatory nodes are (1) the extracellular phosphate concentration x ρ 1 , (2) the concentration of calcitriol, i.e. the active form of vitamin D, x ρ 2 , (3) the PTH concentration x ρ 3 and (4) the FGF23 concentration x ρ 4 .Calcitriol, PTH and FGF23 are hormones that regulate bone metabolism, kidney reabsorption and intestinal absorption of calcium and phosphate, respectively (Melmed et al. 2015, Ch 29).
Following the abstract formulation introduced in Sect.2.3, we can describe this dynamical system by the following system of ODEs Fig. 5 The network associated to the dynamical system (3.1) and its associated specialized subnetworks.a The core network associated to (3.1).From a modeling perspective, I 1 represents the calcium intake, I 2 the phosphate intake, x ι 1 the intestinal concentration of free calcium, x ι 2 the intestinal concentration of free phosphate, x o the extracellular concentration of calcium, x ρ 1 the extracellular phosphate concentration, x ρ 2 the concentration of calcitriol, x ρ 3 the PTH concentration and x ρ 4 the FGF23 concentration.b The I 1specialized subnetwork G I 1 .The I 1 -absolutely super-simple nodes are in green and the appendage nodes in orange.c The I 2 -specialized subnetwork G I 2 .The I 2 -absolutely super-simple nodes are in green and the other I 2 -absolutely simple nodes in pink.Since G I 2 has no appendage nodes, it does not support appendage homeostasis.Consequently, the system does not support pleiotropic-appendage homeostasis.Moreover, since the absolutely super-simple nodes of G I 1 and G I 2 are different from each other (with the exception of o), by Theorem 4.8, the network does not support pleiotropic-structural homeostasis.Consequently, the network supports only coincidental homeostasis.All the nodes belong to both I M -specialized subnetworks G I M , for M = 1 and M = 2 Since we our aim will be to describe the possible types of homeostasis supported by this system rather then to analyze the specific values that the dynamical system can assume, the abstract formulation of the system as given in (3.1) will be enough to our purposes.We refer the reader interested in the precise formulation of the system to Granjon et al. (2017).The multiple-input single-output network associated to the dynamical system above is given by Fig. 5.
The homeostasis determinants with respect to the inputs I 1 and I 2 , respectively, are given by det where the blocks B I 1 and B I 2 are given by and Hence, the vector determinant is given by Note that by the structure of h in (3.3), it is clear that the coordinates of h have no common factor.Hence, according to Definition 2.9, the system does not support pleiotropic homeostasis.As by Proposition 2.8, the system must support in general some type of infinitesimal homeostasis, we conclude that it may present coincidental homeostasis.By our algorithm described in Sect.2.8, this conclusion can also be derived directly from the analysis of the network corresponding to the dynamical system (3.1)(see Fig. 5).
We shall now list all the possible types of coincidental homeostasis that may happen.To simplify notation, we list the factor that appear in the coordinates of h that may be equal to 0, and the corresponding classification of homeostasis.
The theoretical results above give the list of all possible homeostasis types of the general admissible system (3.1).As often, in a model-independent approach, we can not say much about what happens in a specific model equation of the form (3.1), such as the original model in Granjon et al. (2017).However, it may happen that some of the homeostasis types above do not occur in a specific model equation.For instance, it is easy to check if case (1) above can occur in a specific model equation: it is enough to compute f o,x 1 and f ρ 1 ,x ι 2 and verify that they never vanish.When this is the case, we can conclude that this homeostasis type cannot occur in that specific model equation.On the other hand, if they both can vanish, then it may be possible to find homeostasis points by numerical computation.
In general, there is no obstruction for a 'generic' admissible system to display all possible homeostasis types.Moreover, it is not difficult to numerically find a point of infinitesimal homeostasis in a 'generic' admissible system.In Fig. 6 we present the result of a numeric computation of the input-output function x o (I 1 , I 2 ) of a generic admissible vector field (3.1) truncated up to quadratic order.The numeric computation allows us to find that infinitesimal homeostasis occurs at (I 0 1 , I 0 2 ) ≈ (2.9, 12.7).The plateau is located at x o (2.9, 12.7) ≈ 0.09.Near the singularity the function x o (I 1 , I 2 ) is topologically equivalent to a hyperbolic saddle-a Morse singularity in R 2 with normal form h(I 1 , I 2 ) = I 2 1 − I 2 2 (see Golubitsky and Stewart (2018) for more details).Recall that a Morse singularity has codimension 0, thus it is structurally stable, namely any small perturbation of x o (I 1 , I 2 ) (induced by a small perturbation of the admissible vector field) is topologically equivalent to the unperturbed function.The flatness of the input-output function x o (I 1 , I 2 ) near the homeostasis point is reflected in the graph Fig. 6a, which shows that for (I 1 , I 2 ) ∈ [6, 15] × [−2, 7], the value of x o stays in [0, 2].

Reduction to the core network
In this section, unless explicitly stated, we assume that G is a multiple inputs network with input nodes ι 1 , . . ., ι n , and inputs I 1 , . . ., I N .
The definition of core subnetwork implies that the admissible system of Eqs.(2.17) can be written as Now we can freeze the variables x d at an appropriate value and obtain an admissible system for G c from system (4.1).
Lemma 4.1 Suppose that the point X ) is a linearly stable equilibrium of (4.1).Then the G c -admissible system obtained from (4.1) by freezing x d at x * d , given by Proof Clearly, X * c is an equilibrium of (4.2).As shown in Madeira and Antoneli (2022, Lem 3.1), it is linearly stable.Indeed, the Jacobian matrix J of (4.1) evaluated at X * is Therefore, the eigenvalues of J are the same eigenvalues of f d,x d , f u,x u and of the matrix J c , where Since J c is the Jacobian matrix of (4.2) calculated at X * c , it follows that if X * is a linearly stable equilibrium then so it is X * c .
Theorem 4.2 Let x o (I) be the input-output function of the admissible system 2.17 for the network G and let x c o (I) be the input-output function of the admissible system (4.2) for the corresponding core network G c .Then x c o exhibits infinitesimal homeostasis at I * if and only if x o exhibits infinitesimal homeostasis at I * .
Proof For each weighted homeostasis matrix H M , we have: Hence, for each 1 ≤ M ≤ N , we have: where From Lemma 4.1, we have Applying (4.6) and (4.8) to (2.13), we get: Therefore, x o and x c o have exactly the same critical points.

Classification of pleiotropic homeostasis types
In this sub-section, unless explicitly stated, we assume that G is a core multiple inputs network with input nodes ι 1 , . . ., ι n , and inputs I 1 , . . ., I N , with n, N ≥ 2.
We shall now study the subnetworks associated to pleiotropic homeostasis blocks.Bearing this in mind, we start by extending the classification of nodes from Madeira and Antoneli (2022).

Pleiotropic-appendage homeostasis
To study the structure of pleiotropic-appendage homeostasis, we shall first generalize the concepts of path equivalence and appendage subnetworks employed in Wang et al. (2021), Madeira and Antoneli (2022) to the current context.Definition 4.2 Let K be a nonempty subnetwork of G.We say that nodes ρ i , ρ j of K are path equivalent in K (or K-path equivalent) if there are paths in K from ρ i to ρ j and from ρ j to ρ i .A K-path component is a path equivalence class in K. ♦ (c) The appendage subnetwork A G is the subnetwork of G composed by nodes which are I M -absolutely appendage, for all M = 1, . . ., N , and the arrows connecting such nodes.That is, Now we can characterize the structure of pleiotropic-appendage homeostasis.Let B be a pleiotropic appendage block.By a similar argument employed in Madeira and Antoneli (2022), we conclude that B must be the jacobian matrix of the corresponding subnetwork K B .

Pleiotropic-structural homeostasis
Now we shall study the pleiotropic-structural blocks.Let V G be the set of nodes of G, V G ι the set of nodes that are ι m -super simple, for all m = 1, . . ., n and V G I the set of nodes that are I M -absolutely super-simple, for all M = 1, . . ., N .In Madeira and Antoneli (2022), we introduced the notion of absolutely super-simple nodes with respect to the input nodes.This suggests that we can define absolutely super-simple nodes with respect to the inputs.This leads to the question: Which subset of V G is more suitable to base the characterization of pleiotropic-structural subnetworks:V G ι or V G I .The simple, yet paramount, observation that the answer to this question is that both sets are equal.

Lemma 4.5 Let G be a multiple inputs core network. Then
Then, there is at least one input node ι m such that ρ is not an ι m -supersimple node.As G is a core network, there is an input I M such that f I M ,ι m ≡ 0, which implies that ρ is not I M -absolutely super-simple and hence ρ ∈ V G \V G I .On the other hand, if ρ / ∈ V G I , then there exists M such that ρ is not I M -absolutely super-simple ⇒ ρ is not ι m -super-simple, for some m such that f The importance of Lemma 4.5 is that it allows us to study the set V G ι = V G I through either the characterization with respect to the input nodes or to the inputs, whichever is more convenient.In particular, we can easily extend many of the results obtained in Madeira and Antoneli (2022).
A slightly modification of the argument of Lemma 4.5 shows that the set of nodes that are ι m -simple, for all m = 1, . . ., n, and the set of nodes I M -absolutely simple, for all M = 1, . . ., N , are also equal.These observations justify the generalization of the concept of absolutely simple and absolutely super-simple nodes.
Similarly to Wang et al. (2021); Madeira and Antoneli (2022), we say that two elements ρ k > ρ k+1 of V G ι are adjacent when ρ k+1 is the first element of V G ι which appears after ρ k in that ordering.We can now use this concept to introduce the elements that are crucial to characterise pleiotropic-structural homeostasis blocks.Definition 4.7 Let ρ k > ρ k+1 be adjacent elements of V G ι .An ι m -absolutely simple node ρ is between ρ k and ρ k+1 if there exists an ι w o-simple path that includes ρ k to ρ to ρ k+1 in that order, for some w = 1, . . ., n. ♦ The idea is to construct the structural subnetworks employing the concepts above, as it was done in Wang et al. (2021); Madeira and Antoneli (2022).
The absolutely super-simple subnetwork, denoted L(ρ k , ρ k+1 ), is the subnetwork whose nodes are absolutely simple nodes between ρ k and ρ k+1 and whose arrows are arrows of G connecting nodes in L(ρ k , ρ k+1 ).As we can characterise the absolutely super-simple and absolutely simple nodes (and consequently the absolutely supersimple subnetworks) with respect to each input node, we can construct the basic unit of pleiotropic-structural homeostasis in the same way the basic unit of structural homeostasis was constructed in Madeira and Antoneli (2022).Definition 4.8 Let ρ k and ρ k+1 be adjacent absolutely super-simple nodes in G.The absolutely super-simple structural subnetwork L (ρ k , ρ k+1 ) is the input-output subnetwork consisting of nodes in L(ρ k , ρ k+1 ) ∪ B, where B consists of all absolutely appendage nodes that are C S m -path equivalent to nodes in L(ρ k , ρ k+1 ) for some ι m o-simple path S m , for some m ∈ {1, . . ., n}.That is, B consists of all A G -path components B i that are C S m -path equivalent to nodes in L(ρ k , ρ k+1 ) for some S m , for some m ∈ {1, . . ., n}.Arrows of L (ρ k , ρ k+1 ) are arrows of G that connect nodes in L (ρ k , ρ k+1 ).Note that ρ k is the input node and that ρ k+1 is the output node of L (ρ k , ρ k+1 ).♦ We shall employ the characterisation of super-simple structural subnetworks with respect to each of the input nodes.This was the strategy employed in Madeira and Antoneli (2022), and hence we will be able to apply directly the results contained in Madeira and Antoneli (2022, Sect. 3.4.2) to the case of networks with multiple inputs.
First, for ρ k and ρ k+1 adjacent ι m -super-simple nodes in the core subnetwork G m , define as in Madeira and Antoneli (2022) the ι m -super-simple structural subnetwork L m (ρ k , ρ k+1 ) as the input-output subnetwork consisting of nodes in L m (ρ k , ρ k+1 ) ∪ B m , where B m consists of all ι m -appendage nodes that are C m S m -path equivalent to nodes in L m (ρ k , ρ k+1 ) for some ι m o-simple path S m .As usual, arrows of L m (ρ k , ρ k+1 ) are arrows of G m that connect nodes in L m (ρ k , ρ k+1 ).
We notice that (Madeira and Antoneli 2022, Lemma 3.21) is still valid in the current context.Hence, we obtain the following.Proof If B is an irreducible pleiotropic-structural block, then it is a structural block associated to each specialized subnetwork G I M .Fix an input I M and consider the corresponding specialized subnetwork G I M .By Madeira and Antoneli (2022, Thm 3.22), (Golubitsky and Wang 2020, Thm 6.11) and Lemma 4.7, this implies that there exist ι m -absolutely super-simple nodes ρ k M and ρ k M +1 such that K B = L m (ρ k M , ρ k M +1 ) for all m such that ι m is an input node of the specialized subnetwork G I M .Now, as the input and output nodes of all these networks must be the same, we conclude that there exist absolutely super-simple nodes ρ k , ρ k+1 such that for all m = 1, . . ., n, K B = L m (ρ k , ρ k+1 ).By Lemma 4.7, this means that K B = L (ρ k , ρ k+1 ).
The argument in the proof of Theorem 4.8 suggests that, as in the case of inputoutput networks with only one input (Madeira and Antoneli 2022;Golubitsky and Wang 2020), a multiple inputs input-output network supports pleiotropic-structural homeostasis whenever there are more than one absolutely super-simple node.Theorem 4.9 If G has absolutely super-simple nodes other than the output node, then each absolutely super-simple structural subnetwork corresponds to a pleiotropicstructural homeostasis subnetwork.

Conclusion and outlook
In this paper, we present a framework for the analysis and classification of homeostasis types multiple input single output networks.We accomplish this by generalizing and extending the results of Wang et al. (2021) and Madeira and Antoneli (2022) for the classification of homeostasis types in single-input networks single-output networks.Wang et al. (2021) treat the case where the single input parameter affects a single input node and Madeira and Antoneli (2022) consider the case where the single input parameter may affect multiple input nodes.
In the terminology of Golubitsky and Stewart (2022) our theory is an example of a model independent approach.This means that the classification results obtained here provide a complete list of possible behaviors, with respect to homeostasis, that is independent of the model equations-the list depends only on the topology of the network.Which of those behaviors will be observed in a particular realization of the dynamics (e.g. a model equation) depends on the specific form of the dynamics.
We illustrate the application of the theory in several examples.In Sect.2.9, Example 2.15, we analyze the simplest class of multiple inputs networks: the two inputs, three node networks, where each input node is affected by exactly one input parameter (see Fig. 3).In Example 2.16 we have a two inputs, three node network violating this condition-namely, with more than one input parameter affecting the same input node (see Fig. 4).Finally, in Sect. 3 we consider a biologically realistic model for the co-regulation of calcium and phosphate (Granjon et al. 2017) (see Fig. 5).
In three of the eight networks in Fig. 3-cases (b), (e) and (f)-and the network in Fig. 4, we found the 'simultaneous occurrence' of infinitesimal homeostasis and steady-state bifurcation, for certain coincidental homeostasis types (see Remark 2.17).Strictly speaking, this kind of behavior is forbidden by definition, because at a bifurcation point the input-output function becomes ill-defined.However, it is possible, in certain situations, to extend the definition of input-output function to allow for the presence of singular points (see Duncan et al. 2018;Duncan and Golubitsky 2019).These extensions of the notion of homeostasis open up the door for a rich variety of phenomena.For instance, in Mulukutla et al. (2014) the authors investigate glycolysis metabolism and discover a switch mechanism based on a bistability phenomena occurring simultaneously with homeostasis.
The systematic blending of homeostasis and steady-state bifurcations seems to be a promising research avenue.In this regard, the infinitesimal homeostasis approach has some benefit due to its singularity theoretic flavor and the fact that there exists a mature theory of bifurcations based on singularity theory (Golubitsky and Schaeffer 1985;Golubitsky et al. 1988).In fact, Duncan and Golubitsky (2019) is, in part, an attempt to explain the observations of Mulukutla et al. (2014) using singularity theory to uncover the 'interaction' between homeostasis and steady-state bifurcations.
In our examples it seems that the interaction between homeostasis and steady-state bifurcations is 'caused' by the overlapping of the subnetworks associated to certain coincidental blocks in distinct components of the vector determinant.This is distinct from the phenomena discovered in Duncan et al. (2023), where it is shown that an interaction between homeostasis and steady-state bifurcations may occur already in single input node, single input parameter networks.the classification of homeostasis types, we were able to completely characterize the pleiotropic homeostasis types and have provided necessary and sufficient conditions for its occurrence (Sect.4).The main result essentially says that pleiotropic homeostasis types are exactly the homeostasis types that occur in single input parameter, single input node networks.
As for the coincidental homeostasis type, the situation is much more complex.On one hand, were able to obtain some sufficient conditions for its occurrence (see Proposition 2.9).On the other hand, we have given examples where only pleotropic types can occur (Proposition 2.6) and examples where only coincidental types can occur (Example 2.16, cases (a), (b), (c), (d), (f) and the network for calcium and phosphate homeostasis).Furthermore, by Proposition 2.8 any core network must have at least one homeostasis type.Which implies that if there is no coincidental homeostasis type the all homeostasis types must be pleiotropic.All these considerations suggest that a necessary and sufficient condition for occurrence of coincidental types seems rather elusive (see Remark 2.14) and is an important open problem at the moment.

Fig. 1
Fig.1The possible connections in G. Here, inputs are highlighted by dotted circles of distinct colors.The arrows from each input to distinct input nodes are of the same color (the same color of the corresponding dotted circle)

Figure 1
Figure1shows the types of connections which can be found in G. Definition 2.5 Let G be a multiple inputs network.The core subnetwork G c of G is the subnetwork whose nodes are: (i) the input nodes ι 1 , . . ., ι n , (ii) the regulatory nodes σ that are upstream from the output node and downstream of at least one input node, and (iii) the output node o.The arrows of G c are the arrows of G connecting the nodes of G c .♦ .20) Moreover, det(J D I M ) = 0 over J and the irreducible factors of det(J D I M ) never are irreducible scalar factors of h.Proof In case D I M = ∅, the result follows from the convention that det(J D I M ) ≡ 1.In case D I M = ∅, the vestigial subnetwork is composed by nodes that are not downstream from the input nodes affected by the parameter I M .Hence, we can apply to G I M the 'core network' theorem for networks with multiple input nodes and a single input parameter(Madeira and Antoneli 2022, Thm 3.2).The statement about the irreducible factors of det(J D I M ) follows from an argument similar to the one employed inMadeira and Antoneli (2022, Prop 3.8).

I
(a) Input counterweight if B contains partial derivatives with respect to inputs (the simplest counterweight block is of the formf ι m ,I M ), (b) Appendage if B has self-couplings, (c) Structural if B has exactly − 1 self-couplings.♦It follows from the argument in(Madeira and Antoneli 2022, Sec.3.4), applied to the specialized subnetworks G I M , that each homeostasis block of G is of one of the classes in Definition 2.11.Definition 2.12 Let G be a core multiple inputs network.(a) We say that pleiotropic homeostasis occurs when at least one pleiotropic block has vanishing determinant at some fixed input value.The pleiotropic blocks determine the pleiotropic homeostasis types of G. (b) We say that coincidental homeostasis occurs when a N -tuple of coincidental blocks N has simultaneously vanishing determinants at some fixed input value.The N -tuples of coincidental blocks determine the coincidental homeostasis types of G. ♦ Definition 2.13 Let G be a core multiple inputs network and B be a homeostasis block.The homeostasis subnetwork K B associated to B is defined as follows.The nodes of K B are the nodes σ and ρ of G such that f σ,x ρ is a non-zero entry of B. The arrows of K B are the arrows σ → ρ of G such that σ, ρ ∈ K B with σ = ρ.♦

Fig. 3
Fig. 3 Three-node core networks with one input per input node

Fig. 4
Fig. 4 Three-node core network with input affecting both input nodes

Fig. 6
Fig. 6 Numerical computation of the input-output function x o (I 1 , I 2 ) of a generic admissible vector field (3.1) truncated up to quadratic order.Infinitesimal homeostasis occurs at (I 0 1 , I 0 2 ) ≈ (2.9, 12.7), with x o (2.9, 12.7) ≈ 0.09.Panel a shows the 3D plot of the graph of x o (I 1 , I 2 ).Here, the scale of the z-axis (x o ) is different from the the scale of the other two axes.Panel b shows the contour plot (level curves) of x o (I 1 , I 2 ).Near the homeostasis point the function x o (I 1 , I 2 ) is topologically equivalent to hyperbolic saddle.The input-output function was numerically computed using xppaut (Ermentrout 2002) and plotted using r (R Core Team 2023)

Definition 4 . 1
Let G be a multiparameter core network.(a) A directed path connecting nodes ρ and τ is called a simple path if it visits each node on the path at most once.(b) An ι m o-simple path is a simple path connecting the input node ι m to the output node o.(c) A node is ι m -simple if it lies on an ι m o-simple path.(d) A node is ι m -appendage if it is downstream from ι m and it is not an ι m -simple node.(e) A node is I M -absolutely simple if it is an ι m -simple node, for every m such that f ι m ,I M ≡ 0. (f) A node is I M -absolutely appendage if it is an ι m -appendage node, for every m such that f ι m ,I M ≡ 0. (g) An ι m -super-simple node is an ι m -simple node that lies on every ι m o-simple path.(h) An I M -absolutely super-simple node is a node that lies on every ι m o-simple path, for every m such that f ι m ,I M ≡ 0. ♦ It is immediate that the output node o is an I M -absolutely super-simple node, for all M = 1, . . ., N .

Definition 4. 3
The G-complementary subnetwork of an ι m o-simple path S is the subnetwork C S consisting of all nodes of G not on S and all arrows in G connecting those nodes.♦ Definition 4.4 Let G be a multiparameter core network.(a) For every m = 1, . . ., n, we define the ι m -appendage subnetwork A G m as the subnetwork of G composed by all ι m -appendage nodes and all arrows in G connecting ι m -appendage nodes.(b) For every M = 1, . . ., N , we define the I M -appendage subnetwork A G I M as the subnetwork of G composed by all I M -absolutely appendage nodes and all arrows in G connecting I M -absolutely appendage nodes.That is,

Theorem 4 . 3
Let K B be a subnetwork of G associated with a pleiotropic-appendage block B. Then the following statements are valid: (i) Each node in K B is an I M -absolutely appendage node, for all M = 1, . . ., N .(ii) For every ι m o-simple path S, nodes in K B are not C S-path equivalent to any node in C S\K B , for all m = 1, . . ., n; (iii) K B is a path component of A G .Proof Statements (a) and (b) follow by applying (Madeira and Antoneli 2022, Thm 3.11) to each of the core subnetworks G I M .Statement (c) is proved along the same line as (Madeira and Antoneli 2022, Thm 3.11c).Now we shall verify that the conditions listed in Theorem 4.3 are also sufficient to guarantee the existence of a pleiotropic-appendage homeostasis block.Theorem 4.4 Suppose K is a subnetwork of G such that: (i) K is an A G -path component; (ii) For every ι m o-simple path S, nodes in K are not C S-path equivalent to any node in C S\K j , for all m = 1, . . ., n.Then det(J K ) is an irreducible factor of h.Proof Apply (Madeira and Antoneli 2022, Thm 3.13) to each of the specialized subnetworks G I m .The validity of condition (b) of (Madeira and Antoneli 2022, Thm 3.13) for each specialized subnetwork G I M follows directly of condition (b) of this theorem.It is then enough to prove that K is a path component of A G I M , for all m = 1, . . ., n.As K j is a path component of A G = 1≤M≤N A G I M , then for each M = 1, . . ., N , there is a A G I M -path component T M such that K ⊆ T M .By condition (b), it follows that K = T M , for each M = 1, . . ., N .

Definition 4 . 6
Let G be a multiparamete core network.Define a relation on V G ι = V G I as follows: for any pair of nodes σ, τ ∈ V G ι = V G I , σ = ρ, we write σ > ρ when ρ is downstream from σ by all ι m o-simple paths, for any m = 1, . . ., n. ♦ Lemma 4.6 The relation on V G ι = V G I given in Definition 4.6 is a total order.Proof This result is analogous to (Madeira and Antoneli 2022, Lem 3.15) Consider now the ordered elements of
: (R n ) k → R n determines the influence of the joint state (x 1 , . . ., x k ) on the i-th node state x i .The function H i may depend not only on two node states, but may involve multiple nodes concurrently.Network dynamical systems of the form (2.6) have been considered in the groupoid formalism I M associated to the input I M as an H I M -input counterweight, an H I M -structural or an H I M -appendage block.As the derivatives f ι j ,I M would appear in the expression of det(B) if it was a H I M -input counterweight block, we conclude that det(B) must be either an H I M -structural or an H I Madeira and Antoneli (2022)multiple inputs core network and B be a pleiotropic block of G.(i) Then B is either appendage or structural.(ii)Moreprecisely,B is an appendage (respect.structural)blockifandonlyifit isH I M -appendage (respect.H I M -structural) block, for all M = 1, ..., N .Proof (i) From the results ofWang et al. (2021),Madeira and Antoneli (2022), we see that, for each M = 1, . . ., N , B can be classified with respect to the specialized subnetwork G M -appendage block, for all M = 1, . . ., N .
f ι,I N where det(H ) is the homeostasis determinant of the network G, as a polynomial function of the partial derivatives f j,x i .In particular, Condition (2.8) implies that infinitesimal homeostasis occurs if and only if det(H ) = 0. Proof This is a consequence of Lemma 2.1 and of Eq. (2.16) for det H M when G has a single input node.
b) A node ρ is called absolutely simple if and only if it is an ι m −simple node, for all m = 1, . . ., n. Equivalently, ρ is called absolutely simple if and only if it is an I M -absolutely simple node, for all M = 1, . . ., N .♦ (a) A node ρ is called absolutely super-simple if and only if it is an ι m -super simple node, for all m = 1, . .., n.Equivalently, ρ is called absolutely super-simple if and only if it is an I M -absolutely super-simple node, for all M = 1, . .., N .(