Paper deep dive
Topology Inference for Immune System Networks by Using Cell Amount Data
Yushan Li, Rikard Forlin, Dimos V. Dimarogonas, Petter Brodin
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/10/2026, 1:59:02 AM
Summary
This paper proposes a novel topology inference method for immune system networks using cell amount data from cell-depletion experiments. Addressing challenges like limited data availability and lack of standard analytical models, the authors characterize three key properties of immune cell interactions: state non-negativity, ratio-based convergence, and triple signs of topology weights (promoting, inhibiting, or no influence). They construct a new nonlinear model that satisfies these properties and propose a constrained quadratic programming method to infer the network topology from sparse data pairs. Validation on experimental data demonstrates the method's effectiveness.
Entities (10)
Relation Signals (8)
Yushan Li → affiliatedwith → KTH Royal Institute of Technology
confidence 99% · Yushan Li ... Department of Decision and Control Systems, KTH Royal Institute of Technology, Sweden
Rikard Forlin → affiliatedwith → Karolinska Institutet
confidence 99% · Rikard Forlin ... Department of Women’s and Children’s Health, Karolinska Institutet, Sweden
Topology Inference → appliedto → Immune System Networks
confidence 95% · This paper focuses on inferring the topology of a group of immune cells, based on the collected data from cell-depletion based experiments.
Constrained Quadratic Programming → usedfor → Topology Inference
confidence 94% · Finally, based on the constructed model, we propose a constrained quadratic programming method to infer the topology from limited number of data pairs.
Cell-Depletion Experiments → providesdatafor → Topology Inference
confidence 93% · This paper focuses on inferring the topology of a group of immune cells, based on the collected data from cell-depletion based experiments.
Nonlinear Model → enforcesproperty → Ratio-based Convergence
confidence 90% · Then, we construct a new model with simple structure and analytical convenience, and obtain sufficient conditions for the model to accommodate all three properties.
Nonlinear Model → enforcesproperty → Triple Signs of Topology Weights
confidence 90% · Then, we construct a new model with simple structure and analytical convenience, and obtain sufficient conditions for the model to accommodate all three properties.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Recent years have witnessed the advanced development of topology inference research, which helps elucidate the interaction relationships of components in many biological networks. This paper focuses on inferring the topology of a group of immune cells, based on the collected data from cell-depletion based experiments. The problem is very challenging due to i) the lack of standard analytical models for the cell interactions, and ii) the restrictive data availability determined by the huge experiment and time costs. To address these issues, we first leverage certain common knowledge and observations on the experiments to characterize three properties on the cell amounts during the interaction process: state non-negativity, ratio-based convergence, and triple signs of topology weights. Then, we construct a new model with simple structure and analytical convenience, and obtain sufficient conditions for the model to accommodate all three properties. Finally, based on the constructed model, we propose a constrained quadratic programming method to infer the topology from limited number of data pairs. Validation on experiment data demonstrate the effectiveness of the proposed method.
Tags
Links
- Source: https://arxiv.org/abs/2608.07403v1
- Canonical: https://arxiv.org/abs/2608.07403v1
Trouble viewing inline? Open PDF directly →
Full Text
37,634 characters extracted from source content.
Expand or collapse full text
Topology Inference for Immune System Networks by Using Cell Amount Data Yushan Li Rikard Forlin Dimos V. Dimarogonas Petter Brodin Department of Decision and Control Systems, KTH Royal Institute of Technology, Sweden (email: yushanl, dimos@kth.se). Department of Women’s and Children’s Health, Karolinska Institutet, Sweden (e-mail: rikard.forlin, petter.brodin@ki.se). Abstract Recent years have witnessed the advanced development of topology inference research, which helps elucidate the interaction relationships of components in many biological networks. This paper focuses on inferring the topology of a group of immune cells, based on the collected data from cell-depletion based experiments. The problem is very challenging due to i) the lack of standard analytical models for the cell interactions, and i) the restrictive data availability determined by the huge experiment and time costs. To address these issues, we first leverage certain common knowledge and observations on the experiments to characterize three properties on the cell amounts during the interaction process: state non-negativity, ratio-based convergence, and triple signs of topology weights. Then, we construct a new model with simple structure and analytical convenience, and obtain sufficient conditions for the model to accommodate all three properties. Finally, based on the constructed model, we propose a constrained quadratic programming method to infer the topology from limited number of data pairs. Validation on experiment data demonstrate the effectiveness of the proposed method. keywords: Network systems, topology inference, immune cells, network modeling, consensus. †thanks: This work was supported by the Knut and Alice Wallenberg (KAW) Foundation, and the Swedish Research Council (VR). 1 Introduction Network systems have been widely used to model many biological networks, such as brain neurons, genes, and proteins networks (Barabasi and Oltvai, 2004). Topology inference (or identification) has played an important role to understand the intrinsic interaction relations between different entities in the network. For example, it can be used to reveal the connectome topology on brain neurons (Srivastava et al., 2020), or interpret the regulatory mechanism of cells against cancer (Anastasiadou et al., 2018). In this paper, we focus on inferring the topology of a group of immune cells based on the measured data from cell-depletion experiments. Concerning inferring the topology of network systems from data, numerous works have been developed, e.g., causality-based (Dimovska and Materassi, 2021), vector autoregressive based (Zaman et al., 2020), and graph signal processing based (Leus et al., 2023) methods, to name a few. In recent years, lots of research have made promising progress in the inference of biological processes especially on gene regulatory networks (GRNs) (Badia-i-Mompel et al., 2023). For instance, Tsiantis et al. (2018) proposed an inverse optimal control method to identify the underlying optimality principle from time-series data. Dong et al. (2025) utilized the Sinkhorn’s algorithm to infer the signs of promotion/inhibition relationships in GRNs. Lamoline et al. (2025) designed an optimal transport based method to fit a differential equation model and infer GRNs. Despite of the fruitful advances achieved by these works, there is still much left that is difficult to infer. Specifically, for the immune system network, the difficulties of having both a large biological variation between individuals and technological variation between datasets become evident. This makes it difficult to train current models without a large volume of data, a common feature for the aforementioned works. Furthermore, models that can take a less granular overview of cell-cell dependencies from sparse data are lacking, as many focus on GRNs within each cell. Another obstruction that hinders us from inferring the immune network is that the true interaction mechanism of these cells has not yet been fully elucidated. Therefore, different from engineering systems that can be described by well-documented dynamical models, there are no universal models for immune networks. Luckily, some properties about the interaction process are at least well acknowledged. For example, the non-negativity of the cell amount, the triple signs of the topology weights, and the convergence of the cell amounts to a stable baseline after a perturbation, have been established (Perelson and Weisbuch, 1997; Gunawardena, 2010). How to construct an appropriate model that can accommodate these critical properties of the interaction process of immune cells is of great importance. Motivated by the Ockham’s razor, it is meaningful that one can begin with using the simplest linear time-invariant models to fit the data with certain performance guarantees. For instance, given a non-negative initial state, the standard consensus model (Olfati-Saber et al., 2007) can guarantee the state evolution of all nodes will converge a common state. The scaled consensus model (Roy, 2015) is further proposed to allow for a ratio pattern in the converging state. However, these models require the topology weights to be non-negative. Altafini (2013) investigated the bipartie consensus model specifically considering negative weights, and provided convergence guarantees. Aalto et al. (2022) considered a stochastic linear model for the gene regulatory network and investigated the identifiablity problem from the mean and the covariance of the state distribution. Nevertheless, they still lack non-negative constraints on the actual state level. Based on the above observations, this paper aims to construct a simple model to interpret interaction processes and design a method for inferring the immune system’s topology. The contributions are summarized as follows. First, based on the common prior knowledge and experiment observations, we formally characterize three properties for the interaction process of immune cells, including the non-negativity of states, the convergence to a relative stable state, and the triple signs of a topology weight. Second, we inherit the structure simplicity of traditional consensus models to construct a new nonlinear model with an appropriate physical meaning for the immune network. Specifically, sufficient conditions concerning the topology weights and bounds of state ratios are obtained, which guarantee that the required three properties can be met. Finally, based on the experiment data, we provide a constrained quadratic programming method to infer the topology. Both numerical simulations and experiments verify the effectiveness of the proposed model and method. The remainder of this paper is organized as follows. In Section 2, the system model is constructed along with feasibility conditions, and the corresponding inference method is also provided. Numerical simulations and experiments based on real-life data are conducted in Section 3. Finally, Section 4 concludes the paper. 2 Modeling for The Immune Network Consider that the immune cell network is described by a gragh =,ℰG=\V,E\, where =1,⋯,nV=\1,·s,n\ is the set of different types of cells and ℰE is the set of connection edges among the cells. Specifically, the edge (i,j)∈ℰ(i,j) indicates that cell j have influence on cell i, and wijw_ij is the connection weight for (i,j)(i,j). Then, W=[wij]i,j=1n∈ℝn×nW=[w_ij]_i,j=1^n ^n× n constitutes the topology matrix among the n types of cells. Notations. In this paper, we denote ℝ>0nR^n_>0 (ℝ≥0nR^n_≥ 0) as the set of all n-dimensional real-value vectors that have positive (non-negative) elements. Let 1 and InI_n be the all-one vector and n-dimensional identity matrix, respectively. The superscript (⋅)⊺(·) denotes the transpose of a matrix, vec(⋅)vec(·) is the vectorized form of a matrix column by column, and ⊗ represents the Kronecker product. Given a matrix equation M1XM2=M3M_1XM_2=M_3 where the matrices are with compatible dimensions, the vectorization of the equation satisfies (M2⊺⊗M1)vec(X)=vec(M1XM2)=vec(M3)(M_2 M_1)vec(X)=vec(M_1XM_2)=vec(M_3). 2.1 Principles of Experiments and Data Acquisition In the conducted experiments on investigating the dependencies of immune cells, we need to first knock-out a targeted cell type, and then measure the remaining amounts of all cells at some instants. However, due to the huge experiment cost and long time process, we can only collect very few samples for each experimental condition. Specifically, we sample the data at 22 and 2020 hours, respectively, at unstimulated conditions. Notice that the sampled data at each timepoint contain massive information about the interactions, e.g., the population of immune proteins, mRNA and other materials. This paper focuses on the amounts of the cells and uses them to infer the topology. 2.2 Feasibility Analysis of Existing Models Let xix_i be the amount of cell i. Based on the common knowledge and experiment evidences on the immune cell network (Perelson and Weisbuch, 1997; Gunawardena, 2010), we observe the following three properties that the network exhibits on the amount level. • P1): Non-negativity of the state. During the whole interaction process among the cells, the amounts of all cells should be positive, i.e., xi(k)≥0,∀i∈,k≥0. x_i(k)≥ 0,~∀ i ,~k≥ 0. (1) • P2): Convergence to relative stable state. It is acknowledged that for a well-functioned immune network, the cell amounts should remain stable (denoted by x⋆x ) after reacting to counter a virus. Specifically, the amounts of different cells are never the same and thus they have a stable percentage. Mathematically, this property can be formulated as limt→∞x(t)=x⋆(‖x⋆‖2<∞),xj⋆xi⋆=μjμi,∀i,j, _t→∞x(t)=x ~(\|x \|_2<∞),~~ x_j x_i = _j _i,~∀ i,j, (2) where μ∈ℝ>0nμ ^n_>0 is the composition (or relative ratio) profile, satisfying ∑i=1nμi=1 _i=1^n _i=1. • P3): Triple-attribute of the topology weight. Based on clinic research, the influence of one type of immune cell on the other can be roughly classified into three kinds: prompting, inhibition, or none. Mapping these attributes onto the topology weight, it can be formulated as wij>0,if celljprompts celliwij=0,if celljhas no influence on celliwij<0,if celljinhibits celli. \!\! \ aligned &w_ij>0,~&&if cell~j~prompts cell~i\\ &w_ij=0,~&&if cell~j~has no influence on cell~i\\ &w_ij<0,~&&if cell~j~inhibits cell~i aligned .. (3) As discussed in Section 1, the critical limitation in existing linear consensus models lies in the conflict between the state’s non-negativity and the triple signs of a topology weight. To overcome this dilemma, we construct a new model that slightly breaks the model linearity but preserves the listed three properties, which will be analyzed by using nonlinear Perron-Frobenius theory (Lemmens and Nussbaum, 2012). 2.3 The Proposed Model Based on the above arguments, we model the system as the following form x(k+1)=μ⊺x(k)μ⊺Wx(k)Wx(k). x(k+1)= μ x(k)μ Wx(k)Wx(k). (4) Note that this model is not intended as a first-principles mechanistic description of immune regulation. The interaction among immune cells is an extremely complex process that has not been fully understood so far, and it is not our ambition to cover all factors in the immune process. Instead, (4) is a coarse-grained model tailored to the experimental cell-composition data and to the inference objective of this work. The discrete-time index k represents consecutive observation windows in the experiment. Compared with the classic linear mapping x′(k+1)=Wx′(k)x (k+1)=Wx (k) that could represent a fully decentralized interaction, the model (4) further introduces a scaling operation on the state (scaled by μ⊺x(k)/μ⊺Wx(k)μ x(k)/μ Wx(k)) and exhibits certain centralization nature. We observe that this point is reasonable because the immune cell system is commonly regarded to be globally regulated in the human body (Poon and Farber, 2020). Notice that in this model, μ⊺x(k+1)=μ⊺μ⊺x(0)μ⊺Wx(k)Wx(k)=μ⊺x(0), μ x(k+1)=μ μ x(0)μ Wx(k)Wx(k)=μ x(0), (5) which indicates the weighted sum of x is invariant in the iteration process. This invariance property resembles the weighted state sum of a linear consensus process. To ease analysis, we introduce y(k)=x(k)/μ⊺x(k)y(k)=x(k)/μ x(k), and then the model (4) is equivalently written as y(k+1)=F(y(k))=Wy(k)μ⊺Wy(k). y(k+1)=F(y(k))= Wy(k)μ Wy(k). (6) In the subsequent contents, we will mainly focus on model (6) and analyze its convergence. Assumption 1 The state y(k)y(k) is lower bounded by a universal vector β∈ℝ>0nβ ^n_>0, i.e., y(k)≥βy(k)≥β component-wisely. The implication of Assumption 1 lies in two aspects. On the one hand, it indicates that the ratio of a type of cell in the weighted sum of all cells is lowered bounded (i.e., x(k)μ⊺x(k)≥β x(k)μ x(k)≥β). This point is reasonable because in a healthy immune system, a cell type will have a individual-specific lower bound stemming from both inherited and non-inherited effects, otherwise the immune system will not function well. On the other hand, it corresponds to the fact that the cell amounts are always nonnegative. However, since the topology W contains both negative and non-negative entries, it is possible that not all W∈ℝn×nW ^n× n will satisfy Assumption 1. Next, we will demonstrate under what conditions Assumption 1 can be met. Based on the bound vector β, we have μ⊺y(0)=1≥μ⊺βμ y(0)=1≥μ β. Define an auxiliary residual variable r=1−μ⊺β≥0. r=1-μ β≥ 0. (7) Then, we define the restricted μ-simplex set as Δ(β)=y∈ℝ≥0n:μ⊺y=1,y≥β. (β)=\y ^n_≥ 0:~μ y=1,~y≥β\. (8) The following result shows how to guarantee that the mapping F(y)F(y) is invariant under Δ(β) (β). Theorem 1 (Invariance for F) Suppose there exist positive constants blb_l and bub_u such that the following bounds hold for each row of W cl(i)=∑j=1nWijβj+rmin1≤j≤nWijμj≥blcu(i)=∑j=1nWijβj+rmax1≤j≤nWijμj≤bu,blbu≥β. \ aligned c_l(i)&= _j=1^nW_ij _j+r _1≤ j≤ n W_ij _j≥ b_l~\\ c_u(i)&= _j=1^nW_ij _j+r _1≤ j≤ n W_ij _j≤ b_u,~ b_lb_u1≥β aligned .. (9) Then, for all y∈Δ(β)y∈ (β), it holds that bl≤Wy≤bu,F(y)∈Δ(β). b_l1≤ Wy≤ b_u1, F(y)∈ (β). (10) pf First, we prove the state positivity in the dynamic process (6). Since y∈Δ(β)y∈ (β), we decompose y=β+ηy=β+η, where η≥0η≥ 0 by construction. Then, for the i-th element of WyWy, we have (Wy)i (Wy)_i =∑j=1nWij(βj+ηj)=∑j=1nWijβj+∑j=1nWijμj(μjηj). = _j=1^nW_ij( _j+ _j)= _j=1^nW_ij _j+ _j=1^n W_ij _j( _j _j). (11) Notice that r=1−μ⊺β=μ⊺(y−β)=μ⊺ηr=1-μ β=μ (y-β)=μ η, and thus (Wy)i(Wy)_i is bounded by (Wy)i≥∑j=1nWijβj+min1≤j≤nWijμj∑j=1n(μjηj)≥∑j=1nWijβj+rmin1≤j≤nWijμj≥bl, aligned (Wy)_i&≥ _j=1^nW_ij _j+ _1≤ j≤ n W_ij _j _j=1^n( _j _j)\\ &≥ _j=1^nW_ij _j+r _1≤ j≤ n W_ij _j≥ b_l, aligned (12) (Wy)i≤∑j=1nWijβj+max1≤j≤nWijμj∑j=1n(μjηj)≤∑j=1nWijβj+rmax1≤j≤nWijμj≤bu. aligned (Wy)_i&≤ _j=1^nW_ij _j+ _1≤ j≤ n W_ij _j _j=1^n( _j _j)\\ &≤ _j=1^nW_ij _j+r _1≤ j≤ n W_ij _j≤ b_u. aligned (13) Then, we have bl≤Wy≤bub_l1≤ Wy≤ b_u1 and F(y)=Wyμ⊺Wy≥blμ⊺(bu)=blbu≥β, F(y)= Wyμ Wy≥ b_l1μ (b_u1)= b_lb_u1≥β, (14) where the property μ⊺=1μ 1=1 is applied in the second equality. By induction, it follows that y(k)=F(y(k−1))∈Δ(β)y(k)=F(y(k-1))∈ (β) for all k≥0k≥ 0. The proof is completed. □ Theorem 1 gives a sufficient construction for W to ensure that the state is always contained in the cone set Δ(β) (β). Intuitively, (9) has no direct dependence on the real-time state, and indicates that the change from yi(k)y_i(k) to yi(k+1)y_i(k+1) is bounded by (bu−bl)(b_u-b_l). This point corresponds to our common sense that the immune cells amounts will vary in a gradual way assuming a reasonable time-period (Brodin and Davis, 2017), e.g., 2h-20h or even a couple of weeks in between. We then present the following result. Theorem 2 (Convergence of F) Under the conditions of Theorem 1, y(k+1)=F(y(k))y(k+1)=F(y(k)) will converge to a unique fixed point y⋆∈Δ(β)y ∈ (β) satisfying Wy⋆=sy⋆withs=μ⊺Wy⋆. Wy =sy ~with~s=μ Wy . (15) pf To analyze the convergence of the model, we need to borrow some notions from nonlinear Perron-Frobenius theory (Lemmens and Nussbaum, 2012). First, let =ℝ>0nK=R^n_>0 denote the interior of the closed positive cone ℝ≥0nR^n_≥ 0, and define the Hilbert projective metric on K as111This metric is originally defined based on partially ordered vector spaces. Since this paper only focuses on the positive orthant ℝ>0nR^n_>0, we directly give its reduced form here. dH(x,y)=log(maxi,jxiyjyixj),x,y∈. d_H(x,y)= ( _i,j x_iy_jy_ix_j ),~x,y . (16) For a linear operator L satisfying L()⊂L(K) , its projective diameter is defined by δ(L)=supx,y∈dH(Lx,Ly). δ(L)= _x,y d_H(Lx,Ly). (17) Next, we introduce the following Birkhoff’s contraction lemma (Lemmens and Nussbaum, 2014, Theorem 2.9). Lemma 1 If L is a cone-linear mapping with L()⊂L(K) , then L is a contraction in the Hilbert metric, satisfying dH(Lx,Ly)≤κ(L)⋅dH(x,y), d_H(Lx,Ly)≤κ(L)· d_H(x,y), (18) where κ(L)=tanh(δ(L)4)κ(L)= ( δ(L)4 ) is the contraction ratio. Note that the above Lemma was originally targeted at a linear mapping L, and we need to demonstrate how the constructed model y(k+1)=F(y(k))=Wy(k)/(μ⊺Wy(k))y(k+1)=F(y(k))=Wy(k)/(μ Wy(k)) can sufficiently meet the conclusion in Lemma 1. First, notice that i) Δ(β)⊂ (β) is a compact subset in K, and i) Wy⊂Δ(β)Wy⊂ (β) always holds by assumption. Hence, the conclusion (18) directly applies to the linear mapping W, i.e., dH(Wx,Wy)≤κ(W)⋅dH(x,y),∀x,y∈Δ(β). d_H(Wx,Wy)≤κ(W)· d_H(x,y),~∀ x,y∈ (β). (19) Second, as the nonlinear mapping F(y)F(y) only applies a normalization on WyWy, it follows from the definition of dHd_H that ∀x,y∈Δ(β)∀ x,y∈ (β), dH(F(x),F(y)) d_H(F(x),F(y)) =dH(Wxμ⊺Wx,Wyμ⊺Wy) =d_H ( Wxμ Wx, Wyμ Wy ) =dH(Wx,Wy)≤κ(W)dH(x,y), =d_H(Wx,Wy)≤κ(W)d_H(x,y), (20) which means that the normalization does not change the projective direction of the mapping W. Clearly, the mapping F(y)F(y) inherits the same contraction ratio as (19). Finally, since WyWy is constrained by bl≤Wy≤bub_l1≤ Wy≤ b_u1, the projective diameter δ(W)δ(W) is explicitly given by δ(W) δ(W) =supx,y∈dH(Wx,Wy) = _x,y d_H(Wx,Wy) =supx,y∈log(maxi,j(Wx)i(Wy)j(Wy)i(Wx)j) = _x,y ( _i,j (Wx)_i(Wy)_j(Wy)_i(Wx)_j ) =log(bubl)2=2log(bubl)<∞. = ( b_ub_l )^2=2 ( b_ub_l )<∞. (21) Thus, we have κ(W)=tanh(δ(W)4)=tanh(12logbubl)<1. κ(W)= ( δ(W)4 )= ( 12 b_ub_l )<1. (22) By referring to the well-known Banach fixed-point theorem (Latif, 2013), (Δ(β),dH)( (β),d_H) is a non-empty complete metric space with a contraction mapping F(⋅)F(·), and thus F(y(k))F(y(k)) will converge to a unique fixed point y⋆∈Δ(β)y ∈ (β) such that F(y⋆)=Wy⋆μ⊺Wy⋆=y⋆, F(y )= Wy μ Wy =y , (23) which leads to (15) and completes the proof. □ Theorem 2 reveals that if the matrix W satisfies the condition (9), the state x(k)x(k) will converge to a fixed point that is determined by W. This property of F(⋅)F(·) is slightly different from the consensus model with a row-stochastic topology matrix, because the converging state of the latter is also dependent on the initial state. More importantly, considering that the ratio condition (2) is expected to be met, if we suppose y⋆=αμy =αμ, then it follows that W(αμ)=(μ⊺W(αμ))(αμ)⇒Wμ=αμμ, W(αμ)= (μ W(αμ) )(αμ)~ ~Wμ= _μ, (24) where αμ=αμ⊺Wμ _μ=αμ Wμ. Clearly, here μ is a right eigenvector of W. Substituting (24) into y⋆=F(y⋆)y =F(y ), we have y⋆ y =Wμ⊺Wμ=μ⊺μ, = Wμ Wμ= μ μ, (25) x⋆ x =μ⊺μ(μ⊺x⋆)=(μ⊺x(0))μ⊺μμ, = μ μ(μ x )= (μ x(0))μ μ, (26) where the property μ⊺x(k)=μ⊺x(0)μ x(k)=μ x(0) is applied in the second equality. Clearly, the stable cell amount distribution is directly dependent on x(0)x(0) while exhibiting a ratio pattern. Remark 1 Compared with the models M1)-M3), the proposed model has the following advantages regarding the feasibility conditions. i) There are no connectivity requirements and magnitude constraints on the eigenvalues of W for stability concerns, only requiring μ to be a right eigenvector of W by (24). i) It is the cell percentage in the overall weighted amount sum, instead of the amount itself, that requires to bounded, which is more practical for the immune network setting. 2.4 Method Design Under Limited Data With the modeling for the immune cell established, this part shows how to inversely infer the topology from limited data with measurement noises in an optimization framework. Note that the stable reference profile μ and the bound parameters β,bl,bu\β,b_l,b_u\ are given based on our prior knowledge and pre-experiments on the immune system. Considering the measurement noises, denote the data pair at r-th experiment as x~0r,x~1r\ x_0^r, x_1^r\, normalized into y~0r=x~0rμ⊺x~0r,y~1r=x~1rμ⊺x~1r. y_0^r= x_0^rμ x_0^r,~ y_1^r= x_1^rμ x_1^r. (27) We remark that the selection of β is very conservative in practice, and thus y~0r,y~1r≥β y_0^r, y_1^r≥β generally hold in the experiments. If not, we only need to further project y~0r y_0^r and y~1r y_1^r into Δ(β) (β). Note that the data y~1r y_1^r is desired to approximate F(y~0r)F( y_0^r), or equivalently, (μ⊺Wy~0r)y~1r−Wy~0r≈0(μ W y_0^r) y_1^r-W y_0^r≈ 0, which is the residual error of the objective function in our optimization problem. In a vectorized form, this error contained in the data is given by M~D⋅vec(W)≜[M~1⊺,⋯,M~m⊺]⊺⋅vec(W), M_D·vec(W) [ M_1 ,·s, M_m ] ·vec(W), (28) where M~r=[y~1r((y~0r)⊺⊗μ⊺)−(y~0r)⊺⊗In],r=1,⋯,m M_r= [ y_1^r(( y_0^r) μ )\!-\!( y_0^r) I_n ],~r=1,·s,m. Next, we demonstrate how to make W meet the state positivity on Δ(β) (β). Notice that in the condition (9), the function cl(i)c_l(i) is concave while cu(i)c_u(i) is convex regarding W[i,:]W_[i,:]. Hence, cl(i)≥blc_l(i)≥ b_l and cb(i)≤buc_b(i)≤ b_u are all convex constraints. To facilitate solving the problem in a standard quadratic program (QP), we introduce two auxiliary variables pi,qi\p_i,q_i\ and equivalently write (9) as ∑j=1nWijβj+rpi≥bl,pi≤Wijμj,∀j∑j=1nWijβj+rqi≤bu,qi≥Wijμj,∀j, \ aligned & _j=1^nW_ij _j+rp_i≥ b_l,~p_i≤ W_ij _j,∀ j\\ & _j=1^nW_ij _j+rq_i≤ b_u,~q_i≥ W_ij _j,∀ j aligned ., (29) which adds 2n2n more constraints for each row but will not affect the feasibility. Based on Theorem 1, when (29) and y~0r∈Δ(β) y_0^r∈ (β) hold, the properties bl≤Wy~0r≤bub_l1≤ W y_0^r≤ b_u1 and F(y~0r)∈Δ(β)F( y_0^r)∈ (β) will be satisfied automatically. In addition, recall that y⋆=μ⊺μy = μ μ is the fixed point of F(⋅)F(·) regardless of the data, and substituting it into y⋆=F(y⋆)y =F(y ) yields Wμ=sμμ, Wμ=s_μ, (30) where sμ=μ⊺Wμ⊺μs_μ= μ Wμ μ. Hence, (30) should be also treated as a strict constraint. Finally, based on the above formulation, inferring the topology W from limited noisy data pairs is transformed to solving the following convex QP problem minW,pi,qi,sμ \!\! _W,\p_i,q_i\,s_μ~~ ‖M~D⋅vec(W)‖22+γ‖vec(W)‖1 \| M_D·vec(W)\|_2^2+γ\|vec(W)\|_1 (31a) s.t. Wμ=sμμ, Wμ=s_μ, (31b) ∑j=1nWijβj+rpi≥bl,pi≤Wijμj,∀i,j, _j=1^nW_ij _j+rp_i≥ b_l,~p_i≤ W_ij _j,~∀ i,j, (31c) ∑j=1nWijβj+rqi≤bu,qi≥Wijμj,∀i,j, _j=1^nW_ij _j+rq_i≤ b_u,~q_i≥ W_ij _j,~∀ i,j, (31d) where γ>0γ>0 is the regularization parameter associated with ‖vec(W)‖1\|vec(W)\|_1. Remark 2 Note that introducing the additional L1L_1 norm term has two benefits. On the one hand, it could promote a sparse pattern on W, which corresponds to our common knowledge that one type of immune cell is directly influenced by only a few other cells. On the other hand, if only limited amount of experiment data is available (e.g., when m<nm<n), it is very likely that the data matrix M~D M_D has small rank and renders no unique solution when only ‖M~vec(W)‖22\| Mvec(W)\|_2^2 is considered. Considering the introduced L1L_1 norm term and the randomness of M~D M_D in the objective function, the problem will have a unique solution with high probability (see (Tibshirani, 2013, Section 2) for details). 3 Simulations and Experiments 3.1 Numerical Simulations First, we use an example network with 55 nodes, whose topology matrix is given by W=[0.43680.16900.94130.15390.23820.12880.53130.41370.6771−0.041000.4066−0.041801.49690.11050.54940.27690.7737−0.04071.117400.54720.09400.4353]. W= bmatrix0.4368&0.1690&0.9413&0.1539&0.2382\\ 0.1288&0.5313&0.4137&0.6771&-0.0410\\ 0&0.4066&-0.0418&0&1.4969\\ 0.1105&0.5494&0.2769&0.7737&-0.0407\\ 1.1174&0&0.5472&0.0940&0.4353 bmatrix. (32) The ratio profile vector is μ=[0.22,0.15,0.23,0.14,0.26]⊺μ=[0.22,0.15,0.23,0.14,0.26] , and the bound parameters are given by β=/10β= 1/10, bl=0.01b_l=0.01, and bu=10b_u=10. It can be verified that this W satisfies the conditions in Theorem 1 and has a right eigenvector μ. Since W W is identifiable up a scalar ambiguity, we use a revised mean square error (denoted by E(W,W^)E(W, W)) and the ratio of correctly inferred edge signs (denoted by R(W,W^)R(W, W)) to evaluate the inference performance on the topology, E(W,W^)=minα>0‖W−αW^‖Frob2‖W‖Frob2, E(W, W)= _α>0\|W-α W\|_Frob^2\|W\|_Frob^2, (33) R(W,W^)=1−‖sign(W)−sign(W^)‖0n2. R(W, W)=1- \|sign(W)-sign( W)\|_0n^2. (34) The smaller E(W,W^)E(W, W) and higher R(W,W^)R(W, W) indicate better inference performance. Given the initial state x(0)=[345,75,1200,345,457]⊺x(0)=[345,75,1200,345,457] , the simulated results are given in Fig. 1. The evolution of ratios of different components is plotted in Fig. 1(a), where solid and dash lines correspond to the actual and desired ratios, respectively. It is clear that the ratios of the components will converge to μ as the iteration increases. Then, using different amount of data pairs, the topology matrix is obtained by solving the homogeneous equation Mvec(W)=0Mvec(W)=0. As shown in Fig. 1(b), the corresponding E(W,W^)E(W, W) and R(W,W^)R(W, W) are biased when m<5m<5, while perfect when m≥5m≥ 5. This phenomenon matches our intuition that the topology is identifiable up to a scalar when we have an appropriate number of noise-free data. (a) State ratio evolution of different components. (b) Inference error of W W. Figure 1: Inference performance on a numerical example. (a) Cell distribution at 22h and 2020h in one experiment. (b) Inferred topology matrix illustration. (c) Prediction errors on data pairs. Figure 2: Inference performance on real experiments. 3.2 Validation on Real Immune Cell Experiments In the real immune cell experiments (Forlin et al, Unpublished), we consider 1010 types of immune cells: B, Basophil, CD4T, CD8T, DC, Monocyte, NK, Neutrophil, pDC, plasmaB, and order them from 11 to n. The procedures of these experiments are explained in Sec. 2.1. Here we collect m=10m=10 groups of data pairs, and draw the cell amount distribution of one group in Fig. 2(a). Considering the measurement noises, we obtain the topology by solving the optimization problem (31), and visualize it in Fig. 2(b). Note that the values displayed in the matrix blocks are magnified 1010 times for better reading. Since the ground truth topology of the immune network is unknown in practice, we cannot use the metrics E(W,W^)E(W, W) and R(W,W^)R(W, W) to directly evaluate the inference performance. Instead, we adopt the following relative prediction error on a sample Ep(W^,x0)=‖μ⊺x0μ⊺W^x0W^x0−x1‖/‖x1‖. E_p( W,x_0)= \| μ x_0μ Wx_0 Wx_0-x_1 \|/\|x_1\|. (35) Then, the average amount of all components of x1r(r=1,⋯,10)x^r_1(r=1,·s,10) and its corresponding prediction amount are drawn in Fig. 2(c), along with the relative prediction error curve. It is intuitive to find that most of the relative errors are below 0.30.3. Notice that the average of the 10 relative errors is 0.3270.327, with only two of them being larger than 0.30.3. These two results correspond to the cases where the overall cell amount is small. In this regard, the proposed model and inference method achieve acceptable performance on revealing the interaction topology of the immune network. 4 Conclusions In this paper, we investigated the topology inference problem of a class of immune networks. First, we constructed a new nonlinear model to describe the immune network, enjoying the merits of simple structure and physical interpretations. Then, we derived the sufficient conditions for model to guarantee that the state non-negativity and the ratio-based convergence can be achieved simultaneously. Finally, a constrained QP method was provided to infer the topology matrix from data. Numerical simulations and validation on experiment data demonstrated the effectiveness of the proposed method. Future direction includes investigating the topology identifiablity under weak prior parameter assumptions, giving systematic sensitivity analysis to prior parameters, and providing efficient input design under stimulants for practical immune experiments. References Aalto et al. (2022) Aalto, A., Lamoline, F., and Gonçalves, J. (2022). Linear system identifiability from single-cell data. Systems & Control Letters, 165, 105287. Altafini (2013) Altafini, C. (2013). Consensus problems on networks with antagonistic interactions. IEEE Transactions on Automatic Control, 58(4), 935–946. Anastasiadou et al. (2018) Anastasiadou, E., Jacob, L.S., and Slack, F.J. (2018). Non-coding RNA networks in cancer. Nature Reviews Cancer, 18(1), 5–18. Badia-i-Mompel et al. (2023) Badia-i-Mompel, P., Wessels, L., Müller-Dott, S., Trimbour, R., Ramirez Flores, R.O., Argelaguet, R., and Saez-Rodriguez, J. (2023). Gene regulatory network inference in the era of single-cell multi-omics. Nature Reviews Genetics, 24(11), 739–754. Barabasi and Oltvai (2004) Barabasi, A.L. and Oltvai, Z.N. (2004). Network biology: understanding the cell’s functional organization. Nature Reviews Genetics, 5(2), 101–113. Brodin and Davis (2017) Brodin, P. and Davis, M.M. (2017). Human immune system variation. Nature reviews immunology, 17(1), 21–29. Dimovska and Materassi (2021) Dimovska, M. and Materassi, D. (2021). A control theoretic look at Granger causality: Extending topology reconstruction to networks with direct feedthroughs. IEEE Transactions on Automatic Control, 66(2), 699–713. Dong et al. (2025) Dong, A., Georgiou, T.T., and Tannenbaum, A. (2025). Data assimilation for sign-indefinite priors: A generalization of Sinkhorn’s algorithm. Automatica, 177, 112283. Gunawardena (2010) Gunawardena, J. (2010). Models in systems biology: the parameter problem and the meanings of robustness. Elements of computational systems biology, 19–47. Lamoline et al. (2025) Lamoline, F., Haasler, I., Karlsson, J., Gonçalves, J., and Aalto, A. (2025). Dynamic gene regulatory network inference from single-cell data using optimal transport. Bioinformatics (Oxford, England), 41(8), btaf394. Latif (2013) Latif, A. (2013). Banach contraction principle and its generalizations. In Topics in fixed point theory, 33–64. Springer. Lemmens and Nussbaum (2012) Lemmens, B. and Nussbaum, R. (2012). Nonlinear Perron–Frobenius Theory. Cambridge University Press. Lemmens and Nussbaum (2014) Lemmens, B. and Nussbaum, R. (2014). Birkhoff’s version of hilbert’s metric and its applications in analysis. Handbook of Hilbert geometry, 275–303. Leus et al. (2023) Leus, G., Marques, A.G., Moura, J.M., Ortega, A., and Shuman, D.I. (2023). Graph signal processing: History, development, impact, and outlook. IEEE Signal Processing Magazine, 40(4), 49–60. Olfati-Saber et al. (2007) Olfati-Saber, R., Fax, J.A., and Murray, R.M. (2007). Consensus and cooperation in networked multi-agent systems. Proceedings of the IEEE, 95(1), 215–233. Perelson and Weisbuch (1997) Perelson, A.S. and Weisbuch, G. (1997). Immunology for physicists. Reviews of modern physics, 69(4), 1219. Poon and Farber (2020) Poon, M.M. and Farber, D.L. (2020). The whole body as the system in systems immunology. iScience, 23(9), 101509. Roy (2015) Roy, S. (2015). Scaled consensus. Automatica, 51, 259–262. Srivastava et al. (2020) Srivastava, P., Nozari, E., Kim, J.Z., Ju, H., Zhou, D., Becker, C., Pasqualetti, F., Pappas, G.J., and Bassett, D.S. (2020). Models of communication and control for brain networks: Distinctions, convergence, and future outlook. Network Neuroscience, 4(4), 1122–1159. Tibshirani (2013) Tibshirani, R.J. (2013). The lasso problem and uniqueness. Electronic Journal of Statistics, 7, 1456 – 1490. Tsiantis et al. (2018) Tsiantis, N., Balsa-Canto, E., and Banga, J.R. (2018). Optimality and identification of dynamic models in systems biology: An inverse optimal control framework. Bioinformatics, 34(14), 2433–2440. Zaman et al. (2020) Zaman, B., Ramos, L.M.L., Romero, D., and Beferull-Lozano, B. (2020). Online topology identification from vector autoregressive time series. IEEE Transactions on Signal Processing, 69, 210–225.