HETEROGENEOUS CONSENSUS DYNAMICS
THROUGH REACTION-DIFFUSION MODELS ON
GRAPHS/NETWORKS

Manuel Miranda, Pablo Moreno-Spiegelberg, Ernesto Estrada

► To cite this version:

Manuel Miranda, Pablo Moreno-Spiegelberg, Ernesto Estrada. HETEROGENEOUS CONSENSUS DYNAMICS THROUGH REACTION-DIFFUSION MODELS ON GRAPHS/NETWORKS. 2025. hal-04996182

HAL Id: hal-04996182

https://hal.science/hal-04996182v1

Preprint submitted on 18 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# HETEROGENEOUS CONSENSUS DYNAMICS THROUGH REACTION-DIFFUSION MODELS ON GRAPHS/NETWORKS

MANUEL MIRANDA, PABLO MORENO-SPIEGELBERG, AND ERNESTO ESTRADA

**ABSTRACT.** Consensus protocols (linear and nonlinear) are extensively used for multi-agent systems to perform a wide variety of tasks. Its connections with diffusive models make a natural bridge between engineering and social/natural complex systems. These protocols give rise to conservative processes on the graphs representing the systems, and reach homogeneous (consensus) steady states. Here we develop a model that extends the concept of global consensus (diffusion) to dynamical processes which are nonconservative of the graphs, and reach bounded steady states, which are not necessarily a global consensus. That is, they allow heterogeneous steady states in which different subsets of agents can reach a consensus among them, which may be different from the ones reached by other groups. We prove some mathematical results indicating how the structure of the graphs representing the system determines the final states of the dynamics. Finally, we illustrate the applications of the model by considering a heterogeneous rendezvous of a group of agents, as well as the connections of the model with the Bass model for diffusion of innovation and the Lotka-Volterra equations to model the evolution of groups of species that interact in different ways.

## 1. INTRODUCTION

In many engineering systems a group of agents, in contrast with a single one, are organized to execute a mission or perform a task with better operational capability and efficiency. These multi-agent systems (MAS) perform tasks like formation flight of unmanned air vehicles, clusters of satellites, self-organization, automated highway systems, and congestion control in communication networks, among others. In performing these tasks the agents are coordinated through the so-called linear or nonlinear consensus protocols [28, 16, 1, 11, 22]. The MAS is usually represented by a graph $G = (V, E)$ where the set of vertices $V$ represents the agents and the set of edges $E$ represents their interactions, which are encoded in the adjacency matrix $A = [A_{ij}]_{n \times n}$, where $n = |V|$ is the number of vertices in the graph. Then, the state of the agent $i$ at the time $t$ is described by

$$ (1.1) \qquad \dot{u}_i(t) = \sum_{j=1}^{n} A_{ij}(u_j(t) - u_i(t)), $$

in the linear protocol, where $\dot{u}_i(t) = du_i(t)/dt$, or by

$$ (1.2) \qquad \dot{u}_i(t) = \sum_{j=1}^{n} A_{ij}\phi_{ij}\{u_j(t) - u_i(t)\}, $$

<footnote>Corresponding author: mmiranda@ifisc.uib-csic.es.</footnote>

<footnote>Corresponding author: estrada@ifisc.uib-csic.es.</footnote>
in the nonlinear one, where $\phi_{ij}\{\cdots\}$ is an action function obeying a given set of conditions [28]. The models of synchronization on networks, which are widely and extensively studied in the literature [2, 31], belong to this class of nonlinear consensus protocols.

When the linear consensus protocol is expressed in matrix-vector form (see Preliminaries), the model is easily recognizable as the diffusion equation on the graph. Diffusion is a very ubiquitous phenomenon in many real-world complex systems ranging from molecular to social and ecological ones [19]. Then, it is not unexpected that many MAS were originally inspired by the diffusive processes occurring in biological systems in nature, such as the labour between ant colonies, parade of fish schools, formation of bird groups and cooperative hunting of the herd, foraging and migration patterns, etc. [26].

Traditionally, diffusion processes on networks have been modeled by means of conservative models, meaning that the total number of particles in the system remains constant over time [4]. The cases of linear and nonlinear diffusive models expressed before belongs to the class of conservative processes. However, there are many processes on real-world complex systems which are nonconservative on the graphs representing the skeleton of the systems. They include traffic in urban environments due to the emergence and sinking of cars from/to parking spaces which make their number on the streets of a city is not constant over time. Other examples include the diffusive neurotransmission in neuronal systems, where the so-called volume transmission is clearly nonconservative in the neuronal network, the communication in social media like X or Facebook where an user can post a message which can be read by her followers, but also (if not constrained by the user) by non-followers, as well as the diffusive processes in ecology which include births and death of species, the diffusion of microplastics and other materials in the ocean, among others (see [4] and references therein).

Another common characteristic of the (linear and nonlinear) consensus models is the obvious fact that the final steady state of the system is characterized by a global consensus in the state of the agents of the system. That is, when $t \to \infty$, $u_j(t) - u_i(t) = 0$ for every pair of agents in the system. This consensus state is not always observed in real-world systems, as different groups of people can have different opinions, animals can split into diverse groups to perform a task, etc. In the case of engineering systems it is also possible that states in which there is a heterogeneous consensus, i.e., groups of agents performing actions different from other groups, is desirable. Here, we conduct research towards the development of a model which accounts for both, nonconservativeness and heterogeneity in the steady state. We develop our approach on the basis of the Lerman-Ghosh Laplacian operator on graphs [8, 7], which depends on a parameter $\chi \in \mathbb{R}$. It generates nonconservative diffusion when plugged into the graph diffusion equation. However, it possesses a finite non-zero steady state for only a single value of the parameter $\chi$, limiting its applicability in many contexts [4]. Therefore, we propose to develop a logistic-like model incorporating this operator, ensuring that the system's solutions remain bounded, specifically by the value of 1, are nonconservative and not necessarily homogeneous. We start by defining all the concepts and equations used in this work in a Preliminaries section. Then, we proceed to derive the new model through various approaches continue with several applications of this model to the study of heterogeneous rendezvous of agents in the unit square, studying the diffusion of innovations in a social system adapting the Bass model for it, and by analyzing the connections between the new model and the Lotka-Volterra equations for studying interactions between plant pollinators in an ecosystem.
## 2. PRELIMINARIES

Here we consider simple, undirected, connected graphs $G = (V, E)$ with $\#V = n$ vertices and $\#E = m$ edges. Let $A = A(G)$ be the adjacency matrix of $G$. Then, the degree $k_v = \sum_j A_{vj}$ of the vertex $v$ is the number of vertices adjacent to it. As usual in graph theory we will designate $\delta = \min_i k_i$ and $\Delta = \max_i k_i$. Here upper case letters are used for matrices and lower case for vectors and scalar, which will be clarified in the text their nature. Also when comparing two vectors $a$ and $b$, we will use $a \preceq b$ ($a \succeq b$) to indicate that every entry of $a$ is smaller (larger) or equal to every entry of $b$. Let $K$ be the diagonal matrix of vertex degrees. The standard graph Laplacian is defined as $\mathcal{L} := K - A$. The standard diffusion on a graph is expressed by the following equation:

$$ (2.1) \qquad \dot{u}(t) = \frac{du(t)}{dt} = -\gamma \mathcal{L}u(t), u(0) = u^0, $$

where $u(t)$ is a column vector of the mass concentrations at every vertex of the graph at time $t$, $\gamma$ is a parameter indicating the diffusivity of the system (hereafter taken as $\gamma = 1$), and $u(0) = u^0$ is the initial condition. The solution to the abstract Cauchy problem (ACP) is given by:

$$ (2.2) \qquad u(t) = e^{-t\mathcal{L}}u^0, $$

where $e^M = \sum_{k=0}^{\infty} \frac{M^k}{k!}$ is the matrix exponential of the generic matrix $M$. This process is said to be conservative in the graph meaning that $\sum_{i=1}^n u_i(t) = \sum_{i=1}^n u_i^0$ for any $t$. That is, the sum of mass concentrations remains constant in the graph at any time.

A qualitatively different process is provided by the equation:

$$ (2.3) \qquad \dot{u}(t) = -\mathcal{L}_{\chi}u(t), u(0) = u^0, $$

where

$$ (2.4) \qquad \mathcal{L}_{\chi} = \chi I - A, $$

with $I$ being the identity matrix and $\chi \ge 0$ being a parameter.

Although the solution of this diffusion process is given by: $u(t) = e^{-t\mathcal{L}_{\chi}}u^0$, it differs qualitatively from the standard diffusion by the fact that it is not conservative, having solutions that depend on the parameter $\chi$, which are written as [4]:

$$ (2.5) \qquad \lim_{t \to \infty} x(t) = \begin{cases} (\psi_1^T u^0) \psi_1 e^{t(\lambda_1 - \chi)} = \infty & \text{for } \chi < \lambda_1 \\ \left(\sum_j u_j^0 \psi_1(j)\right) \psi_1 & \text{for } \chi = \lambda_1 \\ (\psi_1^T u^0) \psi_1 e^{-t(\chi - \lambda_1)} = 0 & \text{for } \chi > \lambda_1, \end{cases} $$

where $\lambda_j, \psi_j$ are the eigenpair of $A$, where $\lambda_1 > \lambda_2 \ge \dots \ge \lambda_n$. Then, it is clear that $\sum_{j=1}^n x_j(t > 0) \ne \sum_{j=1}^n x_j^0$, such that the mass is not conserved on the vertices of the graph.
Apart from the specific case where $\chi = \lambda_1$ where the system reaches a steady state, in the other two cases the system goes either to infinite mass concentrations at the vertices ( for $\chi < \lambda_1$) or to a situation where such concentrations vanish ( for $\chi > \lambda_1$).

To avoid the before mentioned scenarios of blowing-up or vanishing concentrations at the vertices of the graph, Estrada [4] proposed to consider the Fisher-Kolmogorov-Petrovsky-Piskunov model [6, 15] on the graph using the Lerman-Ghosh Laplacian:

$$ (2.6) \qquad \dot{u}(t) = -\mathcal{L}_{\chi}u + F(u), $$

where $f$ is a function representing the advective force or reaction process, which counters the diffusion process and although the process continues to be nonconservative it reaches a steady state.

In [4] the reaction term was selected such that a fraction of the concentration increased at the vertex $i$ is removed from the node. Because $u_i(t) < 1$ it was considered that the fraction to be removed is equal to $u_i(t)$ multiplied by the amount in which the concentration has increased, i.e.,

$$ (2.7) \qquad \dot{u}(t) = -\mathcal{L}_{\chi}u(t) + u(t) \odot (-\mathcal{L}_{\chi})u(t), $$

where $\odot$ represents the Hadamard or entrywise product between two vectors, i.e., if $a$ and $b$ are two column vectors of the same length $n$, then $a \odot b = [a_1b_1 \dots a_nb_n]^T$. This model has the property that, if $u_0 = u(0) \in [0, 1]$, then $u(t) \in [0, 1] \forall t \ge 0$. Estrada [4] only studied the specific case when $\chi = 0$, which corresponds to the logistic equation:

$$ (2.8) \qquad \dot{u}_i(t) = (1 - u_i(t)) \sum_j A_{ij} u_j(t). $$

In this dynamics, the solution is always $\lim_{t \to \infty} u_i(t) = 1$ for all $i \in V$, which is the typical solution of a logistic model.

### 3. MODEL DEVELOPMENT

Here we will consider the general model:

$$ (3.1) \qquad \begin{aligned} \dot{u}_i(t) &= (1 - u_i(t)) \sum_j (-\mathcal{L}_{\chi})_{ij} u_j(t) \\ &= (1 - u_i(t)) \left( -\chi u_i(t) + \sum_{(i,j) \in E} u_j(t) \right), \end{aligned} $$

for any $\chi \ge 0$.

On the basis of the results previously observed for the nonconservative Lerman-Ghosh diffusion equation, it could be thought that the only possible stable steady states are: $\lim_{t \to \infty} u_i(t) = 1$ for all $i \in V$ if $\chi < \lambda_1$ and $\lim_{t \to \infty} u_i(t) = 0$ for all $i \in V$ if $\chi > \lambda_1$. However, as revealed by the following examples, the complexities of this dynamic overpass significantly that intuition.

To warm up let us consider the the path graph $P_4$, which is a linear chain of 4 vertices labeled in consecutive order from one of the end vertices. The spectral radius of this graph is $\lambda_1(P_4) = \frac{1 + \sqrt{5}}{2}$. We first picked a value of $\chi < \lambda_1$, i.e., $\chi = 0.5$ and performed the simulation of the dynamics with a random initial condition at the vertices of the graph and obtained the plot illustrated in Fig.
<figure><img src="image_2.png" /><figcaption>FIGURE 3.1. Node concentration over time for the Path Graph $P_4$ following the Lerman-Ghosh logistic model. Each subfigure corresponds to different values of $\chi$: 0.5, 3, 1.25 and $\lambda_1$, respectively.</figcaption></figure>

3.1a where as expected $\lim_{t \to \infty} u_i(t) = 1$ for all $i \in V$. We then picked the value $\chi = 3.0$, which obviously obeys $\chi > \lambda_1$. We perform the simulations with the same initial condition as before and as illustrated in Fig. 3.1b we get $\lim_{t \to \infty} u_i(t) = 0$ for all $i \in V$. So far so good. However, when we picked other values of $1 < \chi \le \lambda_1$ for this graph (the bound will be clear later on) and perform simulations with exactly the same initial condition used before, we get completely unexpected results. For instance, for $\chi = 1.25$ and for $\chi = \lambda_1$ we get that the concentrations at vertices 2 and 3 go to an equilibrium value, which is different (and higher) from the one obtained for vertices 1 and 4 (see Fig. 3.1c and 3.1d). That is, if this were a consensus process, for $\chi = 1.25$ and for $\chi = \lambda_1$, the vertices 2 and 3 reach an agreement or a position, which is different from the one obtained between vertices 1 and 4.
To gain more insights into this phenomenon we performed simulations using 200 values of $\chi$ equidistributed between 0 and 3, for which we run the dynamics of equation 2.7 using an Euler scheme and varying the initial conditions, covering all the initial states between 0 and 1 in steps of size 0.1 for each of the nodes. Apart from the graph $P_4$ we also considered the cycle $C_4$, and the star graph $S_{1,3}$. The results are illustrated in Fig. 3.2, which correspond to the mean of the final states of all nodes after $10^4$ timesteps. We recall that the values of $\lambda_1$ for the cycle and star of four vertices are, respectively: $\lambda_1(C_4) = 2$, $\lambda_1(S_{1,3}) = \sqrt{3}$.

Before proceeding let us clarify what we will understand here by the stability of the solutions of the model studied. An equilibrium is said to be *stable in the sense of Lyapunov* if for every $\varepsilon > 0$, there exists a $\delta > 0$ such that any trajectory starting within a $\delta$-neighborhood of the equilibrium remains within the $\varepsilon$-neighborhood for all future times [18], meaning that moving slightly away from the solution still leaves the system close to the solution. While Lyapunov stability guarantees that trajectories remain close to an equilibrium, *asymptotic stability* requires that these trajectories actually converge to the equilibrium as time goes to infinity [14]. More formally, an equilibrium point is asymptotically stable if it is Lyapunov stable and there exists a $\delta > 0$ such that every solution starting in the $\delta$-neighborhood converges to the equilibrium as $t \to \infty$.

As usual we will consider the eigenvalues of the Jacobian, $J = \frac{\partial f}{\partial x}\bigg|_{x=x_e}$, where $\dot{u} = f(u)$ as the criteria of stability. In particular, if none of the eigenvalues of the Jacobian is zero, the fixed point is said to be *hyperbolic* and nonlinear dynamics resemble those of their linearization–according to the Hartman–Grobman theorem [12]. Therefore, fixed points with all the eigenvalues of their Jacobian negative will be asymptotically stable while fixed points with at least one positive eigenvalue of their Jacobian will be unstable. In cases where some eigenvalues have zero real parts, the linear analysis is inconclusive, and higher-order terms must be considered. In a neighborhood of a hyperbolic fixed point, the phase space can be decomposed into invariant subspaces known as the *stable* and *unstable manifolds*. The stable manifold consists of all trajectories that converge to the fixed point as $t \to \infty$, while the unstable manifold comprises trajectories that diverge from the fixed point in forward time [9]. This can be summarised as, for such points, some directions are stable while the others are not. Here we refer to “direction” in the following way. In a similar way as $\frac{\partial f}{\partial x}$ is the change of $f$ in the direction of $x$, here we refer to $du_i/du_j$ as the change of the vertex $i$ in the direction of $j$.

Let us designate $u_i(\infty) := \lim_{t \to \infty} u_i(t)$. Then, the following are general observations of the behavior of the mean final state of the nonconservative logistic model in $P_4$, $C_4$, and $S_{1,3}$:

(1) For $\chi < \delta$, $u_i(\infty) = 1$ for all $i \in V$;

(2) For $\chi > \lambda_1$, $u_i(\infty) = 0$ for all $i \in V$. However, there are also several unstable branches appearing for, due to the fact that $u_i^0 = 1$ for some vertex $i$, then $u_i(t) = 1$ for any $t > 0$, so the system falls in a steady state that is only unstable in the direction of that specific vertex;

(3) For $(\delta, \lambda_1)$ there is a stable steady state in $(0, 1)$;

(4) For $\chi = \lambda_1$ there is not a discontinuous transition to $u_i(\infty) = 0$ for all $i \in V$, but several values in the interval between them are reachable, depending on the initial conditions. There is a sigmoidal shape around this point, which is a computational artifact emerging by the fact that the system enters in a transient (quasistable) zone, caused by the proximity to a bifurcation as explained by [13].

Our main goal here is to study the stable solutions in the interval $(0, 1)$. However, to gain insights into them we will start by exploring ruling out the trivial solutions which may appear in
<figure><img src="image_3.png" /><figcaption>FIGURE 3.2. Steady states numerically obtained for different graphs with 4 nodes, respectively the cycle $C_4$, the path $P_4$ and the star $S_{1,3}$.</figcaption></figure>

this interval. It is straightforward to realize that the solutions $u(t) = 1$ and $u(t) = 0$ are steady states of the equation 2.7. To study their stability, we linearize the system around these solutions and calculate the eigenvalues of the corresponding Jacobian. First, we calculate the derivative of each component with respect to every direction:

$$ (3.2) \qquad (J)_{ij} = \frac{\partial \dot{u}_i}{\partial u_j} = \begin{cases} -\chi + 2\chi u_i - \sum_{(i,j) \in E} u_j & \text{if } i = j, \\ 1 - u_i & \text{if } (j,i) \in E, \\ 0 & \text{otherwise.} \end{cases} $$

Thus, when substituting for $u(t) = 0$, we obtain that the Jacobian is exactly $J = -\chi I + A = -\mathcal{L}_\chi$. The solution is stable if and only if all the eigenvalues of the Jacobian are negative. Then, because the eigenvalues of $\mathcal{L}_\chi$ are $\lambda_i - \chi$, we conclude that $u(t) = 0$ is stable if $\chi > \lambda_1$. On the other hand, if $u(t) = 1$ then the Jacobian is: $J = K - \chi I$, whose eigenvalues are: $k_i - \chi$. Hence, if $\chi < \delta$, all of eigenvalues are negative, obtaining the range for which 1 is a stable solution. This solution obviously includes the case $\chi = 0$ that was analyzed previously by Estrada [4]. Finally, considering the stability of the two solutions given before and using the fact that: $\delta \le k \le \lambda_1$, with equality if and only if the graph is degree-regular, we conclude that for the interval $(\delta, \lambda_1)$ there is a stable steady state in the interval $(0, 1)$.

Apart from $u(t) = 0$ and $u(t) = 1$, another trivial steady state that should be ruled out is the one corresponding to the eigenvector corresponding to the 0-eigenvalue of $\mathcal{L}_\chi$, if it exists. We have seen that, using the first two solutions ($u(t) = 0$ and $u(t) = 1$) and assuming that the initial condition $0 \le u^0 \le 1$, all values of $u(t)$ must lie in the interval $(0, 1)$. On the other hand, if $\lambda_i$ are the eigenvalues of the adjacency matrix of the graph, then the eigenvalues of $\mathcal{L}_\chi$ are $\lambda_i - \chi$. Thus, if we consider all values $\chi < \lambda_1$ it is true that there is a value of $\chi = \lambda_{j \neq 1}$, which makes $\lambda_j - \chi = 0$. In other words, by considering all possible values of $\chi < \lambda_1$, for any graph there is always an eigenvalue equal to zero. Therefore, any of these zero-eigenvalues are steady states of the equation. On the other hand, they are unacceptable, in the sense that they must have at least one negative entry and then do not obey the condition that $0 \le u(t) \le 1$. Consequently, the only
acceptable solution of this kind is for $\chi = \lambda_1$, and the steady state is given by any multiple of $\psi_1$ (the corresponding eigenvector), whose entries are known as the eigenvector centrality, such that $0 \le c\psi_1 \le 1$ is ensured.

### 4. SEARCHING FOR STABLE SOLUTIONS IN $(0, 1)$

As we have seen, the stable solutions for $\chi \in (\delta, \lambda_1]$ usually are neither 0 nor 1. Thus, we have to further study the stable solutions of the equation 2.7. In other words, we will investigate here how to locate those critical points in the diagrams of “mean final state” vs. “$\chi$”. To be clear, we call here critical points to those values of $\chi$ where there is a change in the solutions of the model. In the first section we ruled out all the trivial steady states, which include the ones dealing with the whole Lerman-Ghosh operator. The rest of these steady states are given by the fact that the term $(1-u)$ is local, which means that $u_i = 1$ is a steady state for a vertex $i$ independently of the values of the rest of the vertices.

The conditions under which $u_i \equiv 1$ can happen are the following ones.

* The initial condition is already $u_i^0 = 1$, which doesn't change with time.

* The only stable steady state for such a vertex and the given value of $\chi$ is 1. As we have seen this is equivalent to set $\chi$ to be smaller than the degree of $i$.

We then propose an algorithm that allow us to find solutions for model under study, available at [23]. For that purpose, let us start by defining the vector $d_i$ to be the number of neighbors of node $i$ for which $u_i(\infty) = 1$. Let $G'$ denote the graph obtained by removing from the original graph $G$ all vertices for which $u_i(\infty) = 1$. Let us denote by $\mathcal{L}_\chi(G')$ the Lerman-Ghosh Laplacian of the modified graph. Then, it can be easily seen that finding the steady states for the remaining nodes is equivalent to solving:

$$ (4.1) \qquad \mathcal{L}_\chi(G)u_i(\infty) = 0 \longleftrightarrow \mathcal{L}_\chi(G')u_i(\infty) = d. $$

It has been previously demonstrated that non-trivial steady states emerge for $\chi$ values within the range between $\delta$ and $\lambda_1$. We have observed by computational experiments the following result which we have not been able to prove:

**Conjecture 1.** Let $u$ and $v$ be two steady states of the equation and let $\chi_u < \chi_v$ be the critical values of $\chi$ for each solution. That is, solution $u$ only exists if $\chi \ge \chi_u$. Let $i$ and $j$ be two vertices on the graph such that $k_i < k_j$. Then, if $u_i(\infty) = 1$ ($v_i < 1$) it implies that $v_i(\infty) = 1$ ($u_i(\infty) < 1$), and $u_i(\infty) = 1$ implies $u_j(\infty) = 1$.

This conjecture shows that stable solutions for higher-degree nodes emerge at later times than those obtained for lower-degree nodes. Thus, based on this we develop the following algorithm for any given graph, which allows to find solutions of the model in an iterative way:

(1) Start by assigning the value $u_i(\infty) = 1$ to nodes with the highest degree, which is similar to forcing them to obey the first condition in 4.

(2) Subsequently, remove these nodes to obtain one or more subgraphs.

(3) For each subgraph, write its Lerman-Ghosh Laplacian and apply (4.1) to determine the solution for $u$.

The value of $\chi$ at which this solution is valid is determined by ensuring that $u(\infty) \le 1$. However, it is sufficient to verify this condition only for the vertex with the highest degree in each subgraph, as it is the last vertex for which the solution satisfies this constraint. This process is then repeated recursively for each subgraph, setting the highest degree nodes to 1 in each step. If the conjecture
<figure><img src="image_4.png" /><figcaption>FIGURE 4.1. (a) Illustration of a simple graph used in several examples in this work. (b) Indication of the first step of the reduction of the graph illustrated in (a) by removing the vertex 1 to give rise to three equivalent subgraphs formed by the vertices {2, 5}, {3, 6} and {4, 7}. (c) The result of removing the vertex 2 in the subgraph {2, 5}. The subgraphs enclosed by dotted line are the one for which the equations need to be solved.</figcaption></figure>

results to be false then the algorithm continues to be a valid method for finding the solutions of the model but they can also include steady states which are not the stable ones. To understand how the algorithm works we give here an example.

**Example 2.** Let us consider the graph illustrated in Fig. 4.1(a). The spectral radius of the adjacency matrix is $\lambda_1 = 2$. Then, for any $\chi < \lambda_1$ the vertex labeled as 1, which has degree equal to 3, will have steady state $u_1 = 1$. To find the solution of the model we need to solve a system of 10 equations. However, assigning the value $u_1(\infty) = 1$ is equivalent to remove that vertex from the graph, such that the graph illustrated in Fig. 4.1(b) emerges. Here, there are three equivalent subgraphs formed by the vertices {2, 5}, {3, 6} and {4, 7}.

Due to the equivalence of the three subgraphs we only need to solve the system for one of them, e.g., for {2, 5}, which is equivalent to solve: $Xu(\infty) = y$ where $X = -L_\chi(G')$ and $y$ is the vector that counts how many of the removed vertices are connected to the vertices that remain in $G'$. In this case, $G'$ is the subgraph formed by {2, 5, 6}, and we have

$$ (4.2) \qquad \begin{pmatrix} -\chi & 1 \\ 1 & -\chi \end{pmatrix} \begin{pmatrix} u_2(\infty) \\ u_5(\infty) \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} $$

such that we can obtain $u(\infty) = X^{-1}y$ giving

$$ (4.3) \qquad \begin{pmatrix} u_2(\infty) \\ u_5(\infty) \end{pmatrix} = \begin{pmatrix} \frac{\chi}{\chi^2 - 1} \\ \frac{1}{\chi^2 - 1} \end{pmatrix}. $$

We now determine for which values of $\chi$ we get $u_2(\infty) < 1$ (because 2 is the vertex of the highest degree among the three remaining ones). That is, we solve $1 = \chi/(\chi^2 - 1)$, which is
<figure><img src="image_5.png" /><figcaption>FIGURE 4.2. Illustration of the change of the final state of the dynamics with the heterogeneous consensus model for the graph in Fig. 4.1(a), where the two zones identified in the example are colored.</figcaption></figure>

$\chi^2 - \chi - 1 = 0$ and give $\chi = \{\pm\varphi\}$, where $\varphi = \frac{1+\sqrt{5}}{2}$. Therefore, for $\chi \ge \varphi$ the solutions are $u_1(\infty) = 1$, $u_2(\infty) = u_3(\infty) = u_4(\infty) = \frac{\chi}{\chi^2 - 1}$ and $u_5(\infty) = \dots = u_{10}(\infty) = \frac{1}{\chi^2 - 1}$.

Now, we remove the vertex of the largest degree of the remaining ones, i.e., vertex 2, and we get the structure in Fig. 4.1(c). This indicates that we have to solve the system for the vertex 5. Indeed we proceed by making $u_2(\infty) = u_3(\infty) = u_4(\infty) = 1$, which gives rise to $0 = (1 - u_5(\infty))(-\chi u_5(\infty) + 1)$ and because we can discard the solution $0 = (1 - u_5(\infty))$ we get $u_5(\infty) = 1/\chi$. This obviously give that the critical value for this solution is $\chi = 1$.

In closing we get:

$$ (4.4) \qquad \begin{array}{ll} \text{Zone I} & \left. \begin{array}{l} u_1(\infty) = u_2(\infty) = u_3(\infty) = u_4(\infty) = 1 \\ u_5(\infty) = \dots = u_7(\infty) = \frac{1}{\chi} \end{array} \right\} \text{for } 1 \le \chi \le \varphi, \end{array} $$

and

$$ (4.5) \qquad \left. \begin{array}{l} u_1(\infty) = 1 \\ u_2(\infty) = u_3(\infty) = u_4(\infty) = \frac{\chi}{\chi^2 - 1} \\ u_5(\infty) = \dots = u_{7(\infty)} = \frac{1}{\chi^2 - 1} \end{array} \right\} \text{for } \varphi \le \chi \le \lambda_1. $$

This is corroborated in Fig. 4.2 where the two zones of stability found previously are shadowed in color.

The primary advantage of this approach is that the graphs are divided into smaller subgraphs as nodes are eliminated, reducing the size of the matrices involved and thus accelerating computations. While this assumption is necessary for implementing an ordered recursive process, it also
simplifies the determination of $\chi$ values by limiting the analysis to the highest-degree nodes at each step instead of evaluating all nodes. Also, based on numerical experiments we have the following conjecture.

**Conjecture 3.** Let $i \in V$ with degree $k_i$. Then, a non-trivial solution to the model 2.7 emerges for $\min\{k_i - 1, \lambda_1\} < \chi \le k_i$, where the first inequality is strict if the minimum value is not $\lambda_1$.

The main challenge in deriving a closed-form expression for any graph lies in the dependency of the steady states on both the topology of the graph and on the initial conditions of the dynamical system. Consequently, obtaining a closed formula for $\chi$ values would require an explicit characterization of $(\mathcal{L}_{\chi}(G'))^{-1}$ for every subgraph $G'$ of the original graph $G$. This is inherently unachievable for arbitrary graphs, due to the impossibility of describing explicitly the determinant of the adjacency matrix for any arbitrary graph. We find however some results that help in task of finding explicit solutions for some particular cases.

**4.1. Explicit solutions for low-degree nodes.** Following a contrary approach to the previous iterative procedure, we can start from the bottom and analytically obtain the solutions of special structures inside the graph that obey again the conditions 4. This is specially useful for the lowest degree vertices in the graph, which will be the first ones at which the solution $u_i = 1$ loses its stability. The problem with this approach comes from the multiple ways in which the lowest degree vertices can be connected between them when there are many of such vertices in the graph.

For that purpose, we define the boundary of a set of nodes $\partial S$ as the set of nodes $j \in G \setminus S$ such that exist at least one $i \in S$ so that $i \equiv j$.

We can start with “isolated” vertices, in the sense that every node in their boundary fulfills the conditions 4. We then have the following result.

**Proposition 4.** Consider a simply connected graph with $n > 2$ and $x$ nodes of degree 1. If all nodes of degree $k_i > 1$ satisfy $u_i = 1$, then the steady state $u_i$ for any degree-1 node satisfies:

$$ (4.6) \qquad u_i(\infty) = \frac{1}{\chi}. $$

This solution is valid if $\chi \ge 1$.

*Proof.* Due to the fact that all the low degree vertices are connected to vertices with higher degree, we can assume that $u_i(\infty) = 1$ for any node $i$ such that $k_i > 1$. Hence, those vertices with $k_i = 1$ satisfy the equation $1 - \chi u_i(\infty) = 0$, so it can be solved as $u_i = \frac{1}{\chi}$. This result is valid as long as $0 \le u_i(\infty) \le 1$, so it is satisfied only if $\chi \ge 1$. $\square$

In Figs. 4.3b and 4.3c, the solutions near $\chi = 1$ correspond to these kinds of solution, with mean equal to $\frac{1+\frac{n-1}{\chi}}{n}$ for the star graphs and mean equal to $\frac{n-2+\frac{2}{\chi}}{n}$ for the path.

This reasoning can be used also in the cases in which the minimum degree is higher than 1 if and only if all those vertices with minimum degree are not connected between them as proved in the following result.

**Proposition 5.** Consider a simply connected graph with minimum degree $\delta \ge 1$, without connections among them. If all nodes of degree $k_i > \delta$ satisfy $u_i(\infty) = 1$, then the steady state $u_j$ for any vertex $j$ having $k_j = \delta$ satisfies:

$$ (4.7) \qquad u_i(\infty) = \frac{\delta}{\chi}. $$

This solution is valid if $\chi \ge \delta$.
<figure><img src="image_6.png" /><figcaption>FIGURE 4.3. Stable steady states for the 4-cycle (a), 4-star (b) and 4-path (c) graphs. Red vertical line corresponds to $\lambda_1$.</figcaption></figure>

This could be generalized even if the nodes are not of minimum degree, but all the nodes in their boundary fulfill 4.

**Proposition 6.** Consider a simply connected graph and let $i$ be a node degree $k_i$, such that all of its nearest neighbors have strictly higher degree. Then, for $\chi \ge k_i$, the solution $u_i(\infty) = \frac{k_i}{\chi}$ is a steady state of the equation 2.7.

*Proof.* Simmilarly to the previous proofs, for $\chi < k_i$, node $i$ and all of its neighbors can only have the value of $u_j(\infty) = 1$ as steady state. Then, considering the steady state equation for node $i$ for $\chi \ge k_i$, we obtain:

$$0 = (1 - u_i(\infty))(-\chi u_i(\infty) + k_i);$$

because the $k_i$ neighbors of vertex $i$ have higher degrees than $i$ itself, thus the solution $u_j(\infty) = 1$ is still stable for them.

Hence, excluding $u_i \equiv 1$, the other possible solution is $u_i(\infty) = \frac{k_i}{\chi}$. Finally, the solution is valid as long as $0 \le u_i(\infty) \le 1$. This implies that $\chi \ge k_i$ again. $\square$

On the other hand, if the vertices of low degree are connected among them, other steady states emerge depending on how many of them are interconnected and on the connectivity pattern of such interconnections. As an example, let us consider a graph with two vertices $i$ and $j$ having $k_i = k_j = \delta = 2$, which are connected between them as well as to other vertices with higher degree. Then, following the same steps as before, we assume the trivial state $u_k = 1$ for all the other vertices and thus the vertices $i$ and $j$ satisfy the equation $0 = 2 + u_j - \chi u_i(\infty)$. Assuming the symmetric solution $u_i(\infty) = u_j(\infty)$, we obtain that the solutions are $u_i(\infty) = \frac{2}{\chi-1}$, which is not the same as the state $\frac{2}{\chi}$ that they would have if they were not connected. Due to the fact that there could be more asymmetric solutions, stability in these systems might be more difficult to deal with.

Nonetheless, we can obtain a closed formula for a special case, when cliques of nodes are connected inside a graph. But, for proving it, we need the following result.
**Lemma 7.** The inverse of the Lerman-Ghosh operator on the complete graph of size n, $\mathcal{L}_\chi(K_n)$ is given by:

$$ (4.8) \qquad \frac{1}{\chi^2 - (n-2)\chi - (n-1)} \begin{pmatrix} \chi - (n-2) & 1 & \cdots & 1 \\ 1 & \chi - (n-2) & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & \chi - (n-2) \end{pmatrix}. $$

*Proof.* First, let's calculate the determinant of $\mathcal{L}_\chi(K_n)$. For that, we will call $\det(\mathcal{L}_\chi(K_n)) = x_n$, while the determinant of the matrix that is equal to $\mathcal{L}_\chi(K_n)$ but with a $-1$ in the entry $(1, 1)$ will be called $y_n$.

Now, using the Laplace formula for the determinant for both $x_n$ and $y_n$, we notice that

$$ (4.9) \qquad \begin{aligned} x_n &= -\chi x_{n-1} + (n-1)y_{n-1}; \\ y_n &= -x_{n-1} + (n-1)y_{n-1}, \end{aligned} $$

where we have used the property that the determinant of a matrix with two rows changed is the opposite of the determinant of the original matrix.

We then get the following solutions, which correspond to integer sequences A055137 and A007318, respectively, in the Online Encyclopedia of Integer Sequences:

$$ (4.10) \qquad \begin{aligned} x_n &= \chi^n - \sum_{k=1}^{n-2} \binom{n-2}{k} k\chi^k; \\ y_n &= -\sum_{k=0}^{n-1} \binom{n-1}{k} \chi^k. \end{aligned} $$

Also, notice that the adjoint of $\mathcal{L}_\chi(K_n)$ is $x_{n-1}I + y_{n-1}A$ for $A$ the adjacency matrix of $K_n$. Thus, $(\mathcal{L}_\chi(C_n))^{-1} = \frac{x_{n-1}}{x_n}I + \frac{y_{n-1}}{x_n}A$. After simplifying the expression, we get the final result. $\square$

This example guides us insights about what to do when higher degree vertices are interconnected, with the only restriction that they are cliques of vertices which can also be connected externally to other vertices fulfilling any of the conditions in 4. We then have the following result.

**Proposition 8.** Let G be a graph and let $K_r$, $r \le n$, be a clique of r nodes for some $r > 0$, be a subgraph of G. For each vertex $i \in K_r$, we write its degree in G as $k_i = d_i + r$, for some $d_i \ge 0$.

Consider the dynamical system 2.7 on G. If for every $j \in \partial K_r$ one of the conditions in 4, then any node $i \in K_r$ admits the steady state

$$ (4.11) \qquad u_i(\infty) = \frac{(\chi - r + 1)d_i + \sum_{k \in K_r} (d_k - r)}{\chi^2 - (r-2)\chi - (r-1)} $$

*Proof.* Under any of the previous two conditions, the steady state for any $j \in \partial K_r$ is $u_j \equiv 1$. Thus, for all nodes $i \in K_n$ and ruling out the trivial solution $u_i = 1$, the steady states satisfy the condition:

$$ u_i(\infty) = \frac{d_i - r + \sum_{k \in C_n \setminus \{i\}} u_k}{\chi} $$
Finding thus the steady states is equivalent to solve the system $Au(\infty) = (\overrightarrow{d_i - r})$, where A is in fact the Lerman-Ghosh operator of the complete graph of size $n \mathcal{L}_\chi(K_r)$.

Using the definition of the inverse matrix 7, we obtain that:

(4.12)

$$ u(\infty) = A^{-1} \begin{pmatrix} d_1 - r \\ \vdots \\ d_n - r \end{pmatrix} = \frac{1}{\det(\mathcal{L}_\chi(K_r))} \begin{pmatrix} \chi - (r-2) & 1 & \dots & 1 \\ 1 & \ddots & \dots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \dots & \chi - (r-2) \end{pmatrix} \begin{pmatrix} d_1 - r \\ \vdots \\ d_n - r \end{pmatrix} $$

Doing the multiplication yields the final result.

□

This last proof is a particular scenario of the iterative process described in 4.1, which means that the process, although computationally demanding, can be used to obtain the precise solution for every of the nodes by calculating the inverse of the Lerman-Ghosh Laplacian for the nodes with concentration different than one.

## 5. HETEROGENEOUS CONSENSUS IN THE GENERAL CONTEXT OF MAS AND COMPLEX SYSTEMS

Here we consider some examples of the potential applications of the heterogeneous consensus dynamical system. We start by considering a case of MAS rendezvous in which a group of agents with initial positions on a unit square are directed towards given final positions. In the “classical” rendezvous task, based on the linear consensus protocol, the MAS must direct towards a single point in the square. here we consider the case in which the positions of groups of agents is different from each other, i.e., we want to separate the MAS into subgroups. The second example is based on the connections existing between the current model and the Bass model for the diffusion of innovations in a social group. Our main goal is to show how the selective use of the parameter of the current model can emulate different realistic scenarios in the adoption of an innovation. Finally we study the relations between the current model and the Lotka-Volterra equations to study the interactions between plant-pollinators in a real-world ecosystem, where again we are interested in showing the possibilities of the model to simulate several realistic scenarios.

### 5.1. Heterogeneous rendezvous.
As we have seen in the Introduction, the consensus protocol is widely used in several types of engineering tasks. In the standard consensus dynamics the steady state of the system consists on the equilibrium point in which all the vertices have the same state independently on the initial condition. Let us consider here a rendezvous consensus in the plane, in which we start from initial position of the vertices of a graph and describe their trajectories as $t \to \infty$. Here we will consider the graph illustrated in Fig. 4.1(a) as an illustration. For all the experiments we will consider the same initial conditions to allow a comparative analysis between the different cases. Thus, we study the models

$$ (5.1) \qquad \dot{u}_i(t) = \sum_j -\mathcal{L}_{ij} u_j(t), $$

with $\mathcal{L}$ being the standard Laplacian, and

$$ (5.2) \qquad \dot{u}_i(t) = (1 - u_i(t)) \sum_j (-\mathcal{L}_\chi)_{ij} u_j(t), $$
<figure><img src="image_7.png" /><figcaption>FIGURE 5.1. Illustration of the trajectories of vertices of the graph illustrated in Fig. 4.1(a) using a linear consensus rendezvous on a unit square.</figcaption></figure>

with $u_i \in \mathbb{R}^2$ in both cases, which represent the $x, y$ coordinates of the vertex $i$ in the plane. The initial condition $u_i^0 \in \mathbb{R}^2$ represents here random initial positions of the vertices in the plane.

The results of the standard rendezvous dynamics are illustrated in Fig. 5.1 where it is observed that the position of all vertices converge to the position which corresponds to the average of the initial positions in $x$ and $y$ coordinates.

**Example 9.** Here we study the non-homogeneous consensus dynamics for the same graph in Fig. 4.1(a). We recall that $\lambda_1 = 2$, $\Delta = 3$ and $\delta = 1$ for the graph under study. In Fig. 5.2 we illustrate the trajectories and final states of the vertices of the graph studied when using the parameter $\chi = 0.3$ (a), $\chi = 1.2$ (b), $\chi = 1.8$ (c) and $\chi = 3.0$ (d). As expected according to the theory previously detailed in this work, when $\chi < \delta = 1$ (a) and when $\chi > \lambda_1 = 2$ (d) we observe global consensus in which all vertices go to the same spatial position, either the point (1, 1) or (0, 0) for (a) and (d), respectively. However, for $\chi = 1.2$, which is in the range $1 \le \chi \le \varphi$, we observe in Fig. 5.2 (b) a dissensus state in which the vertices 1, 2, 3, and 4 go to the position (1, 1) of the unit square and the vertices 5, 6 and 7 go to the position $(\chi^{-1}, \chi^{-1}) \approx (0.833, 0.833)$. Finally, in Fig. 5.2 (c) we observe a perfect separation of the vertices by their degrees in the steady state. That is, here $\chi = 1.8$, which is in the range $\varphi \le \chi \le \lambda_1$, such that the vertex with degree 3 goes to the position (1, 1) of the unit square, the vertices of degree 2 go to the position $\left(\frac{\chi}{\chi^2-1}, \frac{\chi}{\chi^2-1}\right) \approx (0.8036, 0.8036)$ and the vertices of degree 1 go to the position $\left(\frac{1}{\chi^2-1}, \frac{1}{\chi^2-1}\right) \approx (0.4464, 0.4464)$. These states of heterogeneous consensus (cases (b) and (c)) are not allowed with the standard consensus model which always reaches the steady state given by the average of the initial position for all vertices.

**Example 10.** We now consider a simple manipulation that allow us to direct the vertices of a graph to different positions in the unit plane in spite of their degrees or other topological features.
<figure><img src="image_8.png" /><figcaption>FIGURE 5.2. Illustration of the trajectories of vertices of the graph illustrated in Fig. 4.1 (a) using the heterogeneous consensus rendezvous on a unit square for $\chi = 0.3$ (a), $\chi = 1.2$ (b), $\chi = 1.8$ (c) and $\chi = 3.0$ (d).</figcaption></figure>

For instance, let us consider again the graph in Fig. 4.1(a), but in which we want to distribute the vertices in the unit square into 6 groups: {1}, {2, 5}, {3}, {4}, {6}, {7}. For that we simply add an anchor vertex $v_a$ which is connected to every other vertex in the graph through weighted edges. This node will have fixed state $v_a = 1$, so it acts like a control of the system. The weights of these edges will change the effective degree of the corresponding vertices. For instance, let us use a weight of 6 for the edge (a, 1), the effective degree of vertex 1 is now $\tilde{k}_1 = 9$. We continue by adding the following weights: 5 to (a, 7), 3 (a, 6), 2 (a, 4) and 1 (a, 3) as well as to (a, 5), respectively. The pair (a, 2) is not connected, i.e, it has weight of zero. This is equivalent to transform the adjacency matrix to $\begin{pmatrix} A & \alpha \\ \alpha^T & 0 \end{pmatrix}$, where $\alpha^T = [6 \ 0 \ 1 \ 2 \ 1 \ 3 \ 5]$. Therefore, we get $\tilde{k}_1 = 9$, $\tilde{k}_2 = \tilde{k}_5 = 2$,
<figure><img src="image_9.png" /><figcaption>FIGURE 5.3. Illustration of the trajectories of vertices of the graph illustrated in Fig. 4.1 using an anchor vertex as described in the text.</figcaption></figure>

$\tilde{k}_3 = 3$, $\tilde{k}_4 = \tilde{k}_6 = 4$, $\tilde{k}_7 = 6$. Although we have only 5 of groups according to the effective degrees, because $\tilde{k}_4 = \tilde{k}_6 = 4$, these two vertices will end up at different positions of the unit square due to the fact that they are connected to vertices with different effective degree, i.e., vertex 4 is connected to vertices 1 and 7 with effective degrees of 9 and 6, respectively, while vertex 6 is only connected to vertex 3 which has effective degree of 3. Let us then consider the use of the heterogeneous consensus dynamics with $\chi = 6.5$ and with a random initial condition except for vertex a which is set to have $u^0(a) = 1$, such that it never changes its position in the unit square. As can be seen in Fig. 5.3 the steady state of the process indicates the existence of the six groups of vertices separated at different positions of the unit square. This example illustrates how we can manipulate the graph to obtain different final states in an heterogeneous rendezvous process with many potential applications in the sciences and engineering.

5.2. **Connections with the Bass model.** Bass model [3] was firstly proposed in 1969 to account for the number of the adopters F of a certain innovation in a population along time. The differential equation of this model reads as:

$$ (5.3) \qquad \begin{aligned} \dot{F}(t) &= p(1-F(t)) + qF(t)(1-F(t)) \\ &= p - (p-q)F(t) - qF^2(t), \end{aligned} $$

where p is a parameter called innovation term or innovation coefficient, by which new people adopt the innovation independently of time, and q is the imitation term, which stands for the rate at which people adopts the innovation when encountering an already adopter person.

**Proposition 11.** The heterogeneous consensus dynamical model 2.7 can be interpreted as a Bass model 5.3 with parameters $p = \sum_{j \in \mathcal{N}_i} u_j$ and $q = -\chi$
*Proof.* From our previous logistic equation, we can rewrite the equation for each node *i* as:

$$ (5.4) \qquad \begin{aligned} \dot{u}_i(t) &= (1 - u_i) \sum_{j \in \mathcal{N}_i} u_j - \chi(1 - u_i)u_i \\ &= \sum_{j \in \mathcal{N}_i} u_j - \left( \sum_{j \in \mathcal{N}_i} u_j + \chi \right)u_i + \chi u_i^2 \\ &= p - (p - q)u_i - q u_i^2. \end{aligned} $$

□

Although this equation would mean a reinterpretation of such model, this variables may not mean exactly the same as in the original model.

First, the equation would model not the number of adopters in a single population but in many of them, whose inter relations are given by the graph edges. Let us then interpret the Bass model as a special case of epidemic spreading where a susceptible population (S) can get “infected” by a given innovation. Then, the values of *u* represent the probability of infection (adopting the innovation) of different individuals in a social network, following the rules:

$$ (5.5) \qquad \begin{aligned} S &= \frac{p = \sum_{j \in \mathcal{N}_i} u_j}{I}; \\ S + I &= \frac{q = \chi}{2S}, \end{aligned} $$

and preserving the quantity $S + I = 1$.

Second, the term *p* is not a fixed parameter but depends on the probability of adoption of its neighbors $u_j$. Hence, the innovation parameter would be in fact the amount of peers' influence received by node *i*. Third, the term *q* could take negative values. These negative values are the most interesting ones as they mean that the probability of adopting an innovation can be self-inhibitory. Indeed, negative values of *q* represent a negative word-of-mouth (WOM) [10], so, instead of new people adopting the innovation because some already have adopted it, it means that adopters tend to reject the already used innovation when they encounter someone not using that same novelty.

Thus, our results would show that, for $\chi \le k_{min}$, the WOM effect, even if negative, would not prevent the innovation to spread fully in all populations (nodes). Therefore, when $k_{min} < \chi < \lambda_1$, the least connected populations will start showing groups of people not adopting the innovation. The critical point $\chi = \lambda_1$ will be the “amount” of negative WOM at which a variety of different scenarios may appear, but the ranking of population according to their percentage of adopters would be given by the eigenvector centrality. Finally, for $\chi > \lambda_1$, innovations die out if not fully established in any of the populations. Analyzing the evolution over time for these values of $\chi$ we can see that some percentages of adopters may start increasing and then die out, which would reflect this effect of negative WOM, as adopters start increasing by the effect of neighboring populations but we can interpret that innovation is so unacceptable ($\chi$ large enough) that the same adopters start refusing it.

**Example 12.** Let us consider an example based on the results of an experiment known as the Columbia University Drug Study [21]. In this study researchers collected data on the first subscription of the drug “gammanym” by 31 physicians in several medical establishments. Then we use here the network formed 31 vertices representing doctors which are interconnected by 67 edges representing face-to-face friendship relations, which were obtained by asking them to name three doctors whom they considered to be personal friends. The graph in question displays a spectral radius of $\lambda_1 \approx 5.7751$, $\Delta = 10$ and $\delta = 1$. In the real-world experiment 17 doctors adopted the innovation
after a period of 18 months (with the first ones adopting it in 3 months), and 14 physicians never adopted the innovation. Obviously, the standard consensus dynamics will assume the adoption of the innovation of all the members of the network as it has been confirmed in previous studies [5]. Recently, we have performed social experiments showing that not always the totality of individuals adopt an innovation, although advantageous for them, due to the existence of stubborn individuals [24]. Therefore, to emulate a situation in which part of the individuals adopt the innovation and others do not, we consider the parameter $\chi = 5.37$, which is smaller than $\lambda_1$ and larger than $\delta$. As can be seen in Fig. 5.4 (a) there are several individuals that adopted the innovation or which has a high probability of adopting it after a relatively large time. However, others drop continuously their probability of adoption which will finally ends up with no adoption. The current model is able to simulate several types of scenarios and we emulate now the hypothetical situation in which a negative WOM is present in the system. For that purpose we use here the heterogeneous consensus dynamical model with parameter $6.43 = \chi > \lambda_1$, which represents a situation in which an innovation is spread across the network with strong negative WOM. For relatively small $t$, e.g., $t \le 100$ in Fig. 5.4 (b) we observe that the majority of individuals start to get discouraged of adopting the innovation, and so their probability of adopting the new drug drops continuously since the beginning of the process. However, some individuals start to increase the probability of adopting the innovation, reaching a peak, after in which it starts to drop. However, as can be seen in Fig. 5.4 (b) the probability of adoption of all individuals in the network finally dies out as a consequence of the strong negative WOM imposed to the system.

<figure><img src="image_10.png" /><figcaption>FIGURE 5.4. (a) Evolution of the adoption of an innovation in the network of physicians studied here. In this case we selected a parameter $\delta < \chi < \lambda_1$ which allows that some individuals adopt the innovation while others do not. (b) Time evolution of the diffusion of an innovation with large negative WOM. In this case the innovation is finally not adopted by the population. In both cases, the x-axis represents time, the y-axis is the percentage of adopters and the colors represent different individuals.</figcaption></figure>
5.3. **Connections with the Lotka-Volterra model.** Lotka-Volterra model [17, 30] accounts for a variety of differential equations that model the evolution of groups of species that interact in different ways (mutualism, competition, predator/prey, etc). The generalized Lotka-Volterra model in vectorial form is given by:

$$ (5.6) \qquad \dot{x}(t) = x(t) \odot (r + Mx(t)); $$

where each entry of $x$ represents the percentage of an species present in an ecosystem with respect to the full capacity of that species there. The vector $r$ stands for the growth of each of the species on its own, while the matrix $M$ encloses the relations between different species (including the internal competitions between themselves).

Recently, Mooij et al. [25] use this model by fixing $r = 1$ and $M = -(I + \tau A)$, to give the analytical points at which ecosystems with a large number of species may survive. This approach gives precise values to the connectance (the measure of how much the species interact in the ecosystem) that was proposed by May as early as in 1972 [20] to solve this same question.

We now prove that the heterogeneous consensus dynamical model can describe the same kind of information as the model of Mooij et al. [25] after some transformations.

**Proposition 13.** Consider the generalized Lotka-Volterra 5.6 equation with $r = \frac{1}{\chi}k - 1$ and $M = \frac{1}{\chi}A - I$. Then, the heterogeneous consensus dynamical model 2.7 can be rearranged to that same equation.

*Proof.* In the heterogeneous consensus dynamical model instead of accounting for the saturation of the species, i.e., the percentage of living animal with respect to the maximum value that the ecosystem can hold, we account for the emptiness of that species. Let us then rearrange the model 2.7 by changing variables $v_i = 1 - u_i$, hence obtaining for each node the equation:

$$ (5.7) \qquad \dot{v}_i = \chi v_i (1 - v_i) - v_i \sum_{j \in \mathcal{N}_i} (1 - v_j) = -(k_i - \chi)v_i + v_i \sum_{j \in \mathcal{N}_i} v_j - \chi v_i^2. $$

Changing the timescale to $T = \chi t$, we obtain the equation:

$$ (5.8) \qquad \dot{v}_i = \left(1 - \frac{k_i}{\chi}\right) v_i + \left(\frac{1}{\chi}\right) v_i \sum_{j \in \mathcal{N}_i} v_j - v_i^2. $$

On the other hand, we can expand the generalized Lotka-Volterra equation 5.6 term-wise and include the choice of $r$ and $M$, obtaining

$$ (5.9) \qquad \begin{aligned} \dot{x}_i &= r_i x_i + x_i \sum_{j=1}^{n} M_{ij} x_j \\ &= r_i x_i + x_i \sum_{j \in \mathcal{N}_i} M_{ij} x_j + M_{ii} x_i^2 \\ &= \left(1 - \frac{k_i}{\chi}\right) x_i + \left(\frac{1}{\chi}\right) x_i \sum_{j \in \mathcal{N}_i} x_j - x_i^2; \end{aligned} $$

which is the heterogeneous consensus dynamical model. □

**Example 14.** In this case we consider a small random network formed by 10 species having 38 connections among them. The network has $\lambda_1 \approx 7.67$, $\Delta = 9$ and $\delta = 6$. We then consider that in this ecosystem species collaborate with a rate $\frac{1}{\chi}$ with respect to a self-competitive interaction,
which is set to a value of 1, and a growth rate of $1 - \frac{k_i}{\chi}$. In this example we will consider the value of $\chi = 8.34$. In Fig. 5.5 (a) it is plotted the mean percentage of species survival, where the red line corresponds to $\lambda_1$, so all the non-zero values to the left of it represents “isolated” species surviving only due to their positive growth rate, while the green line represents the value of $\Delta$. Thus, the mean percentage of survival observed between the green and red lines is equal to one, and represents the survival and saturation of species with negative growth rate. Note that, while the mutualism decreases as $\chi$ increases, the intrinsic growth increases. Moreover, this growth is negative if $\chi < k_i$, which explains why all species extinguish if $\chi < k_{min}$. There are two interesting conclusions that we can observe from applying our results to this model interpretation. First is that, as all species survive and saturate their limits when $\chi > \lambda_1$, this mutualist interactions allow for species with negative growth (which alone tend to die) to survive in the ecosystem (see Fig.5.5(b)). This always occurs in the range $\chi \in (\lambda_1, k_{max})$, which has positive measure for any network apart from the degree-regular ones. Second, if some species is not present in the ecosystem (initial condition $v_i(0) = 0$), then the whole ecosystem gets affected and it's needed a very high value of $\chi$ (big self-grow and small interactions) for the rest of species to reach their saturation point. This is, when one species is missing in the ecosystem, it is usually better for the individuals to focus on self-grow rather than in their inter-relationships even though they are positive and while there is a term of self-competition hindering the increase of each species.

<figure><img src="image_11.png" /><figcaption>FIGURE 5.5. (a) Mean percentage of survival in the random network of mutualistic interacting species for different values of $\chi$. (b) Percentage of saturation over time of different species mutualistically interacting in the random graph studied in this example. Each color represents one species, where the solid lines represent species with positive growth rate while the dashed ones represent those with negative growth rate.</figcaption></figure>
**Example 15.** We now consider a real-world example consisting of plant-pollinator interactions in native forest sites with and without domestic cattle from which [29] constructed 8 plant-pollinator interaction networks in four grazed and four ungrazed sites in native Nothofagus dombeyi forest in and around Nahuel Huapi National Park, Argentina. Here we use the aggregated data set grouping the interactions among 144 species in the 8 site locations, which is available in the database [27]. The original data represents a bipartite network $G = (V_1, V_2, E)$ in which $V_1$ is formed by plants and $V_2$ is formed by pollinators, such that edges connect pollinators to the plants they pollinate. Therefore, we consider the $V_2, V_2$ projection of the original network. That is, the network contains 144 pollinators and two pollinators are connected if and only if they pollinate at least one plant in common. The largest connected component of this network has 90 vertices and 1828 edges, indicating it is a very dense network. For the selection of the parameters of the model we consider that the interactions between pollinators are mutualistic, meaning that if two species pollinate the same plant this increases the chances of that plant of getting reproduced and increasing the chances of both species to find food. We do not consider the case in which the food is not abundant and the interaction between the pollinators get competitive. Here we consider the following three scenarios: (i) a random initial saturation including all species in the network, (ii) the same random initial saturation as in (i) but considering that the two most connected species are extinct, i.e., $u_k(0) = 0$ for $k$ being the two most connected vertices, and (iii) the same random initial saturation as in (i) but considering that the two least connected pollinators are extinct. We start the simulations with very low population of all the species, i.e., $u_i(0) \ll 1$ for all $i \in V$. In Fig. 5.6 we illustrate the results corresponding to the three cases studied: (i) continuous line, (ii) broken line, and (iii) dotted line. When all species are interacting in the ecosystem (case (i)) we observe the typical discontinuous transition at the value of $\chi = \lambda_1$, where all the species reach their maximum population size, $u_i(\infty) = 1$, for all $i \in V$. However, when the most connected species of the ecosystem is extinct (case (ii)), we observe a continuous transition beyond the point $\chi = \lambda_1$, instead of the abrupt one previously observed. This means that the system needs to increase the value of the parameter $\chi$ to reach the maximum saturation of the species population, in spite of the fact that it has not changed neither $\lambda_1$ nor $\Delta$. More interesting is the case (iii), where not only the smooth transition to saturation is observed as in the case (ii), but here the system needs even larger values of $\chi$ to reach its species saturation. The reason is that the reproduction rate of each species is given by $1 - k_i/\chi$. Then, because this network has $\lambda_1 \approx 51.6838$ we observe that the reproduction rate of the most connected species at $\chi = \lambda_1$ is $1 - 83/51.6338 \approx -0.6059$, while for the least connected one it is $1 - 2/51.6338 \approx 0.9613$. That is, while the population of the most connected species is reproducing with a negative rate at this value of $\chi$, the one of the least connected species is reproducing significantly. That is, they are the ones that grow the most, independently of their interactions with the rest of species in the network. Therefore, the effect of the least connected species is very large when they are extinct due to the fact that they are the ones that mainly push the system towards saturation.

## 6. CONCLUSIONS

The main conclusion of this work is that the heterogeneous consensus dynamical model developed here can emulate realistic situations in MAS as well as complex social and natural systems. These scenarios give rise to nonconservative, heterogeneous but bounded steady states which are captured by the one-parameter model developed here. We have illustrated here the cases of heterogeneous MAS consensus, the split of social groups into different opinion groups when an innovation is being adopted and the evolution of interactions between plant-pollinators in an ecosystem via the
<figure><img src="image_12.png" /><figcaption>FIGURE 5.6. Stable states of the example ecosystem: with random initial conditions (continuous blue line), when the maximally connected specie is extinct since the beginning (broken red line), and when the two less connected species are not present (dotted green line). The vertical red line corresponds to the spectral radius of the adjacency matrix and the continuous vertical green line to the largest degree.</figcaption></figure>

connections of the current model with the Lotka-Volterra equations. Other areas of applications may include the neurotransmitters diffusion in neuronal networks, the diffusion of news through social media, etc.

**Acknowledgements.** The authors thank financial support from Agencia Estatal de Investigación (AEI, MCI, Spain) MCIN/AEI/ 10.13039/501100011033 and Fondo Europeo de Desarrollo Regional (FEDER, UE) under the María de Maeztu Program for units of Excellence in R&D, grant CEX2021-001164-M). E.E. and M.M. thank Agencia Estatal de Investigación (AEI, MCI, Spain) MCIN/AEI/ 10.13039/501100011033 under grant PID2023-149473NB-I00. M.M. thanks financial support PRE2020-092875 by MCIN/AEI/10.13039/501100011033 and by FSE invierte en tu futuro.

REFERENCES

[1] Abdollah Amirkhani and Amir Hossein Barshooi. "Consensus in multi-agent systems: a review". In: *Artificial Intelligence Review* 55.5 (2022), pp. 3897-3935.

[2] Alex Arenas et al. "Synchronization in complex networks". In: *Physics reports* 469.3 (2008), pp. 93-153.

[3] Frank M Bass. "A new product growth for model consumer durables". In: *Management science* 15.5 (1969), pp. 215-227.

[4] Ernesto Estrada. "Conservative vs. non-conservative diffusion towards a target in a networked environment". In: *The Target Problem*. 2023. URL: https://hal.science/hal-04229218.

[5] Ernesto Estrada and Eusebio Vargas-Estrada. "How peer pressure shapes consensus, leadership and innovations in social groups". In: *Scientific reports* 3.1 (2013), p. 2905.
[6] Ronald Aylmer Fisher. "The wave of advance of advantageous genes". In: *Annals of eugenics* 7.4 (1937), pp. 355-369.

[7] Rumi Ghosh and Kristina Lerman. "Parameterized centrality metric for network analysis". In: *Physical Review E—Statistical, Nonlinear, and Soft Matter Physics* 83.6 (2011), p. 066118.

[8] Rumi Ghosh et al. "Non-conservative diffusion and its application to social network analysis". In: *Journal of Complex Networks* 12.1 (2024), cnae006.

[9] John Guckenheimer and Philip Holmes. *Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*. Springer New York, 1983. ISBN: 9781461211402. DOI: 10.1007/978-1-4612-1140-2. URL: http://dx.doi.org/10.1007/978-1-4612-1140-2.

[10] Mariangela Guidolin and Piero Manfredi. "Innovation diffusion processes: Concepts, models, and predictions". In: *Annual Review of Statistics and Its Application* 10.1 (2023), pp. 451-473.

[11] Muhammad Majid Gulzar et al. "Multi-agent cooperative control consensus: A comparative review". In: *Electronics* 7.2 (2018), p. 22.

[12] Philip Hartman. "On the local linearization of differential equations". In: *Proceedings of the American Mathematical Society* 14.4 (Aug. 1963), pp. 568-573. ISSN: 1088-6826. DOI: 10.1090/S0002-9939-1963-0152718-3. URL: http://dx.doi.org/10.1090/S0002-9939-1963-0152718-3.

[13] Alan Hastings et al. "Transient phenomena in ecology". In: *Science* 361.6406 (2018), eaat6412.

[14] Hassan K Khalil. *Nonlinear Systems: Pearson New International Edition*. 3rd ed. London, England: Pearson Education, Nov. 2013.

[15] AN Kolmogorov, IG Petrovsky, and NS Piskunov. "Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problem". In: *Bull. Moscow State Univ. Ser. A: Math. Mech* 1.6 (1937), pp. 1-25.

[16] Yanjiang Li and Chong Tan. "A survey of the consensus for multi-agent systems". In: *Systems Science & Control Engineering* 7.1 (2019), pp. 468-482.

[17] AJ Lotka. *Elements of physical biology*. Williams and Wilkins, 1925.

[18] A. M. Lyapunov. "The general problem of the stability of motion". In: *International Journal of Control* 55.3 (Mar. 1992), pp. 531-534. ISSN: 1366-5820. DOI: 10.1080/00207179208934253. URL: http://dx.doi.org/10.1080/00207179208934253.

[19] Naoki Masuda, Mason A Porter, and Renaud Lambiotte. "Random walks and diffusion on networks". In: *Physics reports* 716 (2017), pp. 1-58.

[20] Robert M May. "Will a large complex system be stable?" In: *Nature* 238.5364 (1972), pp. 413-414.

[21] Herbert Menzel and Elihu Katz. "Social relations and innovation in the medical profession: The epidemiology of a new drug". In: *Public opinion quarterly* 19.4 (1955), pp. 337-352.

[22] Mehran Mesbahi and Magnus Egerstedt. "Graph theoretic methods in multiagent networks". In: (2010).

[23] Manuel Miranda. *LogisticLerman-Ghosh*. URL: https://github.com/ManuM2B/LogisticLermanGhosh/.

[24] Manuel Miranda et al. "Indirect social influence and diffusion of innovations: An experimental approach". In: *PNAS nexus* 3.10 (2024), pgaec409.

[25] MN Mooij et al. "Stable coexistence in indefinitely large systems of competing species". In: *Proceedings A. Vol. 480*. 2299. The Royal Society. 2024, p. 20240290.

[26] BL Partridge. "The chorus-line hypothesis of maneuver in avian flocks". In: *Nature* 309.6 (1984), pp. 344-345.
[27] Tiago P. Peixoto. *The Netzschleuder network catalogue and repository*. Apr. 2023. DOI: 10.5281/zenodo.7839981. URL: https://doi.org/10.5281/zenodo.7839981.

[28] R Olfati Saber and Richard M Murray. "Consensus protocols for networks of dynamic agents". In: *Proceedings of the 2003 American Control Conference, 2003*. Vol. 2. IEEE. 2003, pp. 951–956.

[29] Diego P Vázquez and Daniel Simberloff. "Changes in interaction biodiversity induced by an introduced ungulate". In: *Ecology Letters* 6.12 (2003), pp. 1077–1083.

[30] Vito Volterra. *Variazioni e fluttuazioni del numero d'individui in specie animali conviventi*. Società anonima tipografica” Leonardo da Vinci”, 1926.

[31] Chai Wah Wu. "Synchronization in networks of nonlinear dynamical systems coupled via a directed graph". In: *Nonlinearity* 18.3 (2005), p. 1057.

INSTITUTE FOR CROSS-DISCIPLINARY PHYSICS AND COMPLEX SYSTEMS (IFISC) CSIC-UIB, PALMA DE MALLORCA, SPAIN.