Limits of multifunctionality in tunable networks.

Nature is rife with networks that are functionally optimized to propagate inputs to perform specific tasks. Whether via genetic evolution or dynamic adaptation, many networks create functionality by locally tuning interactions between nodes. Here we explore this behavior in two contexts: strain propagation in mechanical networks and pressure redistribution in flow networks. By adding and removing links, we are able to optimize both types of networks to perform specific functions. We define a single function as a tuned response of a single "target" link when another, predetermined part of the network is activated. Using network structures generated via such optimization, we investigate how many simultaneous functions such networks can be programed to fulfill. We find that both flow and mechanical networks display qualitatively similar phase transitions in the number of targets that can be tuned, along with the same robust finite-size scaling behavior. We discuss how these properties can be understood in the context of constraint-satisfaction problems.

Many naturally occurring and synthetic networks are endowed with specific and efficient functionality. For example, allosteric proteins globally adjust their conformation upon ligand binding to control the activity of a distant active site (1, 2). Venation networks in the leaves of plants are highly optimized for water and nutrient transport (3). In some cases, networks can change their function depending on the needs of the system; vascular networks in animals (4–6), fungi (7), and slime molds (8) can reroute the transport of fluids, enhancing or depleting nutrient levels to support local growth or activity. Modern power grids must precisely distribute electrical energy generated from a limited number of sources to a large number of consumers with widely varying consumption needs at different times (9). All of these networks are optimized to some degree, by evolution via natural selection, dynamic reconfiguration, or human planning.

A key aspect of such functionality is the complexity of a specific task. We define a “function” as an optimized response of a localized component of a network when another predefined, localized component of the system is activated. A “task” is then defined as the collective response of a set of individual functions due to a single input. The number of functions representing a specific task is the task complexity.

In this work we address the limits of complexity for a single task: How many functions composing a single task can be programed into a network? We consider two examples: (i) mechanical networks—in which nodes are connected by central-force harmonic springs—locally flexing in response to an applied strain and (ii) flow (or resistor) networks—in which nodes are connected by linear resistors—locally producing a pressure drop due to an applied pressure at the source. These systems are related; flow networks are mathematically equivalent to mechanical networks embedded in one spatial dimension—but with a nontrivial node topology (10).

Macroscopic properties of mechanical networks, such as their bulk and shear moduli, can be tuned by modifying only a tiny fraction of the springs between nodes (11–13) [in contrast to random removal (14)]. This idea was extended to show that such networks can be tuned to develop allosteric behavior via selective spring removal (15–17). Allostery in these systems is a single-function task in which a randomly selected spring (the target) responds in a specified way to a strain imposed on a separate pair of nodes (the source). Here we study complex tasks in which multiple targets are controlled by a single source. We study the scaling of the maximal complexity of a task with network size by asking how many individual targets can be successfully tuned simultaneously and show that in both flow and mechanical networks, the limit of task complexity is set by a phase transition.

Network Tuning Protocol

Our method for tuning networks follows the general scheme described in our previous work (15) with slight modifications. We start with 2D configurations of soft spheres with periodic boundary conditions created using standard jamming algorithms. We construct networks by placing nodes at the center of each sphere and links (edges) between nodes corresponding to overlapping particles. This ensemble of networks is used for both spring networks, in which edges are unstretched central-force springs, and flow networks, in which edges are resistive conduits. By using the same set of nodes and edges for both systems, we can directly compare results. We chose this ensemble because it is disordered and provides initial networks with properties reminiscent of the corresponding biological systems. Elastic networks with close-range interactions have often been used to model proteins (18), while many natural flow networks have high numbers of closed loops (19) and are highly interconnected (20). However, our networks do not exhibit modular topologies which can appear in some systems (9). While we briefly touch on modularity, an in-depth study of modular networks is outside the scope of this work.

For each network, a pair of source nodes is chosen randomly, along with a set of target edges. Our goal is to tune the extension (or pressure drop) of each target edge, indexed by , in response to an extension (pressure drop) applied to the source nodes by adding and removing edges from the network. We explore two different types of sources: pairs of nodes connected by a randomly chosen edge and pairs of nodes that are each chosen randomly anywhere in the network (see for global compression and shear sources in mechanical networks).

To control the response of the targets, we define the response ratio for each target. Each is in general a collective property of the network; the response of each target is a function of the total network structure. Before tuning the network, we measure the initial extension (pressure drop) to obtain the initial response ratio of each target . We then tune the response ratio of each target so that its relative change compared with the initial state is greater than or equal to a specified positive constant ; that is, we tune each response ratio to satisfy the constraintThus, for mechanical networks we require contracting edges to contract more and expanding edges to expand more. For flow networks, we require the magnitude of the pressure drop to increase without changing the direction of the flow through each target link.

Our optimization scheme involves minimizing a loss function which penalizes deviations from the constraints in Eq. (). Each optimization step consists of either removing a single link or reinserting a previously removed link to modify the network topology in discrete steps. More specifically, at each step we measure the resulting change in the loss function for each single-link removal or reinsertion and remove or reinsert the link to most decrease the loss function.

Fig. 1 depicts examples of both flow and mechanical networks which have been tuned using our prescribed method for the two different types of applied sources. Fig. 1 shows flow and mechanical networks, respectively, tuned to respond to a source applied to a pair of nodes connected by an edge. Fig. 1 shows the same networks, but with a pair of source nodes that are not connected by an edge.

Fig. 1.

Networks tuned to display multifunctional responses. Each network starts with the same initial topology and the same choice of four target edges (corresponding nodes shown in green). A and C are flow networks while B and D are mechanical networks. For each network a source extension (pressure drop) is applied to a pair of source nodes (shown in red). In A and B the pair of source nodes is connected by an edge, while in C and D the source nodes are not connected by an edge. For flow networks, response ratios are tuned to , while for the mechanical networks they are . The edges removed by tuning are shown as thick blue lines. For flow networks, the magnitude of the pressure on each node is indicated by the node size and the sign of the node pressure is represented by the shape. For mechanical networks, the node displacements are shown as black arrows.

Results

For both flow and mechanical networks, we explore the effects of various aspects of the tuning problem, with particular focus on task complexity. Fig. 2 displays typical results for the fraction of networks that can be tuned successfully, , for flow and mechanical networks, respectively. Data are shown for a randomly chosen edge source and randomly chosen target edges with a desired relative change in target response of . System sizes range from to 4,096 nodes. Each value of is calculated by tuning at least 512 independent randomly generated networks. At low , while at large , drops to zero.

Fig. 2.

(A and B) The fraction of satisfied configurations for (A) flow networks and (B) mechanical networks as a function of number of targets for systems of nodes. Results are shown for a pressure or extension applied to a single source edge with a desired response ratio of . Curves are smoothing splines and estimated error bars are shown for binomially distributed data (). (C) Scaling collapse for all for four cases: flow networks with an edge source (red circles) and with a node pair source (blue triangles) and mechanical networks with an edge source (green squares) and with a node pair source (black triangles). In each case, we plot vs. , where at and is the interval in over which . (D and E) The transition points (D) and width of the transition (E) are reasonably described by power laws in with fits for giving exponents 0.67 and 0.65 for flow networks and 0.71 and 0.74 for mechanical networks with an edge and node pair source, respectively. In the same order, the power-law fits for have exponents of 0.71, 0.66, 0.74, and 0.66.

In Fig. 2, we plot the transition curves for all system sizes for the four cases studied on the same axes. Using the smoothing spline interpolations shown in Fig. 2 (), we estimate the number of targets at which . Next, we estimate the width of the transition, , taken as the interval in over which . We attempt to collapse each curve by plotting vs. . We find a similar functional form for all cases, with only a slight difference between flow and mechanical networks near .

Fig. 2 shows that flow networks and mechanical networks have similar power-law behaviors for and . Both the transition location and width scale approximately as with . Because the scaling exponent for is less than 1, the critical fraction of functions that can be tuned simultaneously approaches zero as goes to infinity, even though the number of simultaneously tuned functions diverges with system size. Thus, small networks are relatively more tunable than large ones. In addition, the sublinear scaling of the transition width shows that drops more rapidly with as increases, implying that the crossover becomes sharp as . At the same time, Fig. 3 shows that the average number of links that need to be removed for a successful tuning operation grows approximately linearly with the number of targets. Thus, those networks that can be tuned successfully typically require only removal of a constant fraction of edges. Together, our results suggest with consistent with zero. For , this implies a random first-order transition in the thermodynamic limit, with a discontinuity in and power-law finite-size scaling. Such hybrid transitions are typical of constraint–satisfaction problems.

Fig. 3.

(A and B) Power-law behavior of the average number of removed edges as a function of number of targets for (A) flow networks and (B) mechanical networks for various system sizes . Included networks correspond to those that have been tuned successfully in Fig. 2 with an edge source and desired change in target response of . Error bars indicate the error on the mean. Power laws with an exponent of 1.0 are depicted as black dashed lines for comparison.

Discussion

We framed the problem of the maximum number of target edges that can be tuned successfully in a mechanical or flow network as a type of discrete constraint–satisfaction problem, in which we asked how many inequality constraints can be satisfied simultaneously. This places the tuning of multifunctionality in the context of a variety of other problems including jamming (21), spin glasses (22), the -SAT problem (23), -core percolation (24), and the perceptron (25). Much progress has been made by linking such transitions to the statistical physics of critical phenomena. The hallmark of these systems is the emergence of a SAT–UNSAT transition between regions in parameter space where the constraints can always (or with high probability) be satisfied and regions where the system is frustrated, such that not all constraints can be satisfied simultaneously (25). In mean-field and in some cases in finite dimensions, the SAT–UNSAT transition is a random first-order transition, with a discontinuous jump in the order parameter (the fraction of satisfied configurations ) as in a first-order phase transition, but with power-law scaling as in a second-order transition.

We have demonstrated a SAT–UNSAT transition in the complexity of a single task that can be tuned into disordered mechanical and flow networks. In both cases, the maximum task complexity diverges with a power law that is sublinear in , the number of nodes in the network. The width of the SAT–UNSAT transition (relative to ) vanishes as diverges, showing that the transition is a true phase transition.

Although we find for the four cases displayed in Fig. 2, both and can vary, depending on a variety of factors. These factors include (i) the local or global nature of the source, (ii) the magnitude of desired change in target response , (iii) disorder in the link topology, (iv) initial coordination of the network, and (v) the choice of whether to tune the link tensions (currents) or extensions (pressure drops) (). The values of lie in the range of 0.6–0.8, with the exception of one case of 1.0 for a very large relative change in target response of (). We find that the behavior is not well described by a power law for tuning negative relative changes in target response () and for tuning small changes in current or tension. The former case is still under investigation, while the latter exception has a simple explanation ().

Overall, the divergence of the maximum number of tunable targets with system size and the corresponding vanishing of the transition width (indicating the existence of a phase transition) are very robust observations for positive and sufficiently large relative changes in target responses. We note also that both mechanical networks and flow networks exhibit very similar quantitative behavior despite the fact that flow networks are purely topological, requiring no explicit spatial embedding.

The SAT–UNSAT transition of the task complexity problem introduced here represents a distinct class of discrete constraint–satisfaction transitions due to a complication that arises in the form of the constraints. When tuning a mechanical network, the removal of links can introduce soft modes, making it impossible to uniquely evaluate the network response and subsequently tune a given target. Similarly, in a flow network the tuning process can lead to regions being disconnected from the source, making it impossible to tune any target in that region. To avoid such cases, at each step of the tuning process we are forced to exclude specific link removals (). In both mechanical and flow networks, we find that it becomes more and more likely to introduce a soft mode/disconnected region as the task complexity increases. This makes the problem more difficult to tackle both numerically and analytically compared with previously studied constraint–satisfaction transitions and may lead to differences in the nature of the transition.

For mechanical functions, a perfectly engineered mechanism (e.g., a pair of chopsticks, which creates a large displacement at the tips in response to strain applied where they are held) may perform exactly one function superlatively well, but we have shown that more complex network structures are able to adapt to a number of functions that diverge with the system size. The same argument holds for flow networks: An optimally engineered distribution network is a topological tree, perfectly suited for a specified task but at the same time “rigid,” in the sense that it cannot easily adapt to other tasks. The networks that we have studied are more complex than a pair of chopsticks or a topological tree, and this allows them to be tuned successfully to perform arbitrarily complex tasks.

Our finding that a disordered network topology allows for tunability may have relevance to real biological networks. For example, the development of certain vascular structures within animals is characterized by the initial appearance of a tightly meshed disordered network of veins (the vascular plexus) that is subsequently pruned and tuned to its function (26). The initial disordered network may be a prerequisite of the great variability and versatility seen in natural networks. The tuned mechanical networks serve as simple models for multifunctional allostery in proteins (with a single regulatory site that can control more than one active site, e.g., refs. 27 and 28) or multifunctional metamaterials. Our flow network results give insight into how to control, for example, blood and oxygen distribution in vascular systems or power in an electrical network. Indeed, we find very similar behavior in a flow network with nonplanar topology derived from the United Kingdom (UK) railroad network, which exhibits a high degree of modularity. exhibits a qualitatively similar transition in the number of targets that can be tuned compared to the transitions in the networks studied here ().

Our results raise a number of issues for future investigation. The divergence in the task complexity and vanishing of the transition width with system size are reasonably well approximated by power laws but may deviate for larger system sizes (). The measured exponents appear to depend on many specific properties of the problems studied. This may be due to corrections to scaling or to a more fundamental deviation from power-law scaling. Also, it is not clear what conditions on the network topology are necessary to observe the transition we see. For example, networks with high degrees of modularity may not be able to support tasks spanning multiple neighborhoods. However, our results for the UK railroad network suggest that even in this case, we observe identical qualitative behavior, but with an overall decreased corresponding to the possibility of choosing sources and targets in different neighborhoods (). More generally, it has not been investigated how the results depend on network structure/topology and dimensionality nor how they depend on the tuning algorithm. For instance, the values of and might be higher for simulated annealing, which explores a wider region of solution space than the minimization algorithm studied here.

One further aspect of our results deserves mention: A simple function that controls only a single pair of target nodes can be achieved in an extremely large number of ways. We have shown that a task can be complex with randomly chosen target nodes controlled by a single source. However, if one is interested in controlling only a single target, one can create different paths for its control by choosing any of the other nodes in the system also to be a target. Likewise, one could specify a third node to be controlled as well, etc. That means that there are at least ways of creating that simple function. Because we find , for this lower bound is smaller than the prediction of solutions in the large- limit (16).

Here we studied the limits of the complexity of a single task. It would be interesting to understand how many different tasks can be designed successfully and whether that is controlled by a similar SAT–UNSAT transition. Finally, we note that for the mechanical and flow networks studied here, the behavior is governed by a discrete Laplacian operator (29)—mechanical networks obey force balance on each node and flow networks obey Kirchhoff’s laws. However, many networks, such as gene regulatory networks, metabolic networks, social networks, etc., are nonconservative. Moreover, the problems we have studied are linear in their couplings but ecological networks or neural networks, for example, are typically nonlinear. It is known that even nonconservative and/or nonlinear networks, such as the Hopfield model and jammed packings, can support SAT–UNSAT transitions as well (30, 31). It would interesting to study systematically how conservation constraints and linearity affect the nature of the transition.

Materials and Methods

Linear Response.

Our networks are described by a set of nodes and edges. The response of a flow network to external stimuli is represented by a pressure on each node . Analogously, the response of a -dimensional mechanical network is the -dimensional displacement vector of each node. Each edge linking nodes and is characterized by either a conductance or stiffness, denoted in both cases. For mechanical networks, , where is the stretch modulus per unit length and is the rest length. Initially, we set all stretch moduli identically to one. Similarly, for flow networks we set all conductivities to one. Removing an edge corresponds to setting to zero, whereas reinserting an edge corresponds to setting back to its original value.

To calculate the response of each type of network, we minimize the corresponding functional. In the case of flow networks, we minimize the power loss through the network,where indicates a sum over all edges. For mechanical networks, we minimize the elastic energywhere is a unit vector pointing from node to node in the undeformed configuration. The power loss for a flow network can be mapped to the energy of a mechanical network for by mapping the pressure on each node to a 1D displacement (10). In this case, the unit vectors are scalars with values of either , which drop out when squared; the embedding of the network in space does not matter as is expected for flow networks.

Minimizing Eq. for a flow network in the presence of externally applied boundary currents on each node , we obtain a system of linear equations characterized by a graph Laplacian ,where is an -dimensional vector of node pressures and is an -dimensional vector of external currents on nodes. We define the vector so that the pressure and current on the th node are and , respectively. Similarly for mechanical networks, minimizing Eq. in the presence of externally applied forces, we obtainwhere is a -dimensional vector of node displacements and is a -dimensional vector of external forces on nodes. Again we define the matrix to pick out the displacement and force on the th node, and . The matrix is the matrix of second derivatives known as the dynamical or Hessian matrix and can be interpreted as graph Laplacian where each element is a matrix. We define the Laplacian, denoted , as a generalized version of the standard Laplacian matrix. The case corresponds to the Laplacian of a flow network (or a 1D mechanical network) such that , while for , is a Hessian for a -dimensional mechanical network; i.e., . The th block of the Laplacian iswhere is nonzero only if edge exists.

Consequently, the response of either type of network is calculated by solving the corresponding set of linear equations rewritten aswhere and are the appropriate -dimensional response and source vectors, respectively. To apply a pressure drop or edge extension source, we use a bordered Laplacian formulation.

Bordered Laplacian Formulation.

Calculating the linear response requires solving Eq. . However, there are two complications. The first is that the Laplacian operator is in general not invertible due to the presence of global degrees of freedom. For a periodic network, in dimensions, there are global translational degrees of freedom. Second, we apply edge extension (pressure drop) sources, rather than tension (current) sources. These sources can be applied as constraints on the system. Using a bordered Laplacian formulation, we add a constraint for each global translation and for the source.

First, we define the extension (or pressure drop) of the source aswith source nodes and . The unit vector points from node to and is a scalar in the case of a flow network. The vector is defined to extract the extension of the source from the full vector of node displacements. We specify the desired extension as . Additionally, we define the vectors for corresponding to translations of the entire system uniformly along the th axis. We define the Lagrangianwhere the parameters and are Lagrange multipliers. We include the Lagrange multipliers as additional unknown parameters that must be determined in our calculations. We find solutions by extremizing the Lagrangian with respect to both the displacements and the Lagrange multiplier. We rewrite the Lagrangian in matrix form:The vector is size with elements and is a size matrix with columns . In this context we can further condense notation, writing the Lagrangian aswhere we define the bordered Laplacian as a block matrix of second derivatives of the Lagrangian:We also define the bordered displacement and force vectors and , respectively, each of size as

As a result, the system of equations we must solve is now . The bordered Laplacian is invertible due to the presence of the constraints and solving this equation is straightforward.

Tuning Loss Function.

Framed according to Eq. , the problem of tuning a complex task can be viewed as a constraint–satisfaction problem. The goal is to find a set of stiffnesses (conductivities) that simultaneously satisfy each constraint in Eq. . To study this problem numerically, we recast it as an optimization problem in the style of ref. 25, in which we define an objective function that penalizes deviation of the system’s behavior from the desired multifunctionality. Thus, we introduce the loss functionwhich is a function of the set of all of the spring constants (conductivities) and is composed of a sum over the set of target edges to be tuned. For each target edge we define the residualwhich measures how close each target is to being tuned successfully. The Heaviside function is included so that if , i.e., the response ratio has increased at least by the desired proportion , then the residual does not contribute to the loss function.

Optimization Method.

Our method for tuning a network involves minimizing the loss function in Eq. . In the spirit of refs. 11 and 15, our optimization consists of removing or reinserting previously removed edges from the network one at a time, modifying the network topology in discrete steps. More specifically, we use a greedy algorithm in which we remove or reinsert the edge which minimizes the loss function at each step. This requires a calculation of the new response for each possible move.

Suppose we have a network whose stiffnesses at the current step are for all valid where some might already be zero, having been removed at previous steps. Our goal is to measure the change in response when the stiffness of edge is changed by an amount . We note that the Laplacian can be decomposed aswhere the equilibrium (or incidence) matrix of size defines the mapping of nodes to edges (29, 32) and is a size diagonal matrix of edge stiffnesses such that . We can define a bordered incidence matrix by appending rows of zeros to , giving us a corresponding decomposition of the bordered Laplacian . The change in response iswith the corresponding change in the bordered Laplacian with the vector . We now need to calculate the inverse of the updated bordered Laplacian. This can be done using the Sherman–Morrison formula (33)The change in response is thenThe new response is then used to calculate an updated loss function.

To reduce numerical error and maintain the numerical invertibility of the bordered Laplacian, we define the quantityIf is less than , we do not remove an edge. This quantity can be shown to be the contribution of an edge to the states of self-stress in mechanical systems (12, 34). By ensuring that every removed edge has some contribution to the states of self-stress, then by Maxwell–Calladine counting, we are guaranteed that no zero modes are introduced (35).

We repeatedly add or remove edges until either the loss function is explicitly zero (i.e., all constraints are satisfied) or the relative change in the objective function is less than .