Fibroblast-Specific Protein-Protein Interactions for Myocardial Fibrosis from MetaCore Network

Abstract

Myocardial fibrosis is a major pathologic disorder associated with a multitude of cardiovascular diseases (CVD). The pathogenesis is complex and encompasses multiple molecular pathways. Integration of fibrosis-associated genes into the global MetaCore network of protein-protein interactions (PPI) offers opportunities to identify PPI with functional and therapeutic significance. Here, we report the generation of a fibrosis-focused PPI network and identification of fibroblast-specific arbitrators driving reparative and reactive myocardial fibrosis. In TGF-$\beta$-mediated fibroblast activation, developed network analysis predicts new regulatory mechanisms for fibrosis-associated genes. We introduce an efficient Erdös barrage approach to suppress activation of a number of fibrosis-associated nodes in order to reverse fibrotic cascades. In the network model each protein node is characterized by an Ising up or down spin corresponding to activated or repairing state acting on other nodes being initially in a neutral state. An asynchronous Monte Carlo process describes fibrosis progression determined by a dominant action of linked proteins. Our results suggest that the constructed Ising Network Fibrosis Interaction model offers network insights into fibrosis mechanisms and can complement future experimental efforts to counteract cardiac fibrosis.

1. Introduction

Myocardial fibrosis is a major pathologic disorder associated with a multitude of cardiovascular diseases (CVD) [1]. In the heart, tissue fibrotic remodeling is characterized by abnormal fibroblast activation and excessive extracellular matrix (ECM) protein accumulation [2]. Although a number of factors have been implicated in orchestrating the fibrotic response, tissue fibrosis is dominated by a central mediator: transforming growth factor-$\beta$ (TGF-$\beta$) [3]. Sustained TGF-$\beta$ production leads to a continuous cycle of growth factor signaling and deregulated matrix turnover [4]. TGF-$\beta$ binds a heteromeric receptor complex composed of type I and type II serine/threonine kinase receptors that activate smad-dependent gene transcription [5]. The converging lines of evidence suggest that activation of the Smad 3 cascade plays an essential role in extracellular matrix protein gene expression and regulates fibrosis tissue deposition in fibrotic remodeling of the infarcted heart [6]. Abnormalities in TGF-$\beta$/Smad3 signaling underlie inflammatory diseases and promote tumor emergence. TGF-$\beta$ is also central to immune suppression within the tumor microenvironments, and several studies have reported roles in tumor immune evasion and poor responses to cancer immunotherapy [7]. In addition to activation of Smad-dependent signaling, the TGF-$\beta$ receptor complex can also modulate TGF-$\beta$-activated kinase 1 (TAK1 or MAP3K7) [8].

However, despite intensive research, the biomolecules that orchestrate fibrosis are still poorly understood and as a result, effective strategies for limiting fibrosis are lacking [2,4]. Considering the complex heterogeneity of fibrosis, research strategy on a system-level understanding of the disease using mathematical modeling approaches is a driving force to dissect the complex processes involved in fibrotic disorders. Recently, we have reproduced the classic hallmarks of aberrant cardiac fibroblast activation leading to fibrosis and provided a powerful toolbox for characterization of cardiac fibroblast activation [9]. Although the pathogenesis of fibrotic remodeling has not been well identified, accumulated evidence suggests that multiple genes/proteins and their interactions play important roles in disease scenario [10].

Fibrotic remodeling is a complex process caused by genetic abnormalities that alter protein-protein interactions (PPI) [9,11]. In the heart, PPI interfaces represent a highly promising, although challenging, class of potential targets for therapeutic options. In fibrosis, PPI form signaling nodes and hubs that transmit pathophysiological cues along molecular networks to achieve an integrated biological output, thereby promoting fibrogenesis and fibrosis progression. Thus, pathway perturbation, through disruption of PPI critical for fibrosis, offers a novel and effective strategy for curtailing the transmission of pro-fibrotic signals. Deciphering of fibrosis-specific PPI would uncover new mechanisms of fibrotic signaling for therapeutic interrogation. Thus a mathematical analysis of the global PPI MetaCore network [12] can provide useful deep insights in the understanding of fibrosis progression processes.

In this work we use the PPI MetaCore network [12] to perform a mathematical modeling of fibrosis progression. The results reported in [9] determined the protein pro-fibrotic responses as a feedback on TGF protein stimulation. Thus these results [9] support the known fact that the TGF protein plays an important role in fibrosis tissue [3] and establish proteins with most positive and most negative response in cardiac fibroblasts.

A variety of biological applications of the MetaCore network are described in [13,14]. The general statistical properties of the MetaCore network are presented in [15]. The applications of Google matrix algorithms to the fibrosis PPI based on the MetaCore network are described in [16] using the fibrosis responses obtained in [9]. We note that the Google matrix algorithms [17,18,19] find a variety of useful applications in modern complex networks [20] including World Wide Web, Wikipedia, world trade etc.

In our opinion an important advantage of the MetaCore network is that it presents a global network structure of PPI with about 40,000 proteins and important molecules. This allows to perform a mathematical modeling of fibrosis progression over the whole network. Thus in this work we present such a modeling which assumes that at an initial stage there is a certain number of proteins which activate a fibrosis progression of other proteins via their network links to other proteins. These initial proteins are always marked as red network nodes and assumed to be always red permanently generating fibrosis progression (thus being always red or with fixed Ising spin up with ${\sigma }_i=1$). Among the other proteins, we assume that there is an initial group of proteins which inhibit fibrosis progression being of fixed blue color (thus being always blue or with Ising spin down with ${\sigma }_i=-1$). The remaining other proteins (or network nodes) are assumed to be in neutral state of white color or spin being zero (${\sigma }_i=0$). The fibrosis progression is modeled by the asynchronous Monte Carlo process when at each step a spin of a given node (which is not fixed) is determined as the total spin sign of other nodes linked to it with certain weights of Markov chain transitions. In a certain sense this rule corresponds to a situation when a person takes an average opinion of his network friends linked to him. After many of such Monte Carlo steps the system converges to a steady-state when all network spins obtain a fixed polarization being up (red or fibrosis activated) or down (blue or healthy protected) and some of the nodes may remain at zero neutral spin state (white) being practically disconnected from the red or blue network nodes.

We note that the above Monte Carlo steps describe a process of opinion formation and the global vote of nodes (society members) between red and blue options. For complex social networks [20] there are numerous studies of opinion formation, reviewed in [21]. Various interesting features of voter models and opinion formation on networks had been obtained and described in [21,22,23,24,25,26,27]. Recently it was argued that the opinion formation on the world trade network can be linked with country preference to trade in one or another currency (e.g., US dollar or hypothetical BRICS currency) [28]. The important new element appeared in [28] is that opinion of certain network nodes (countries) is considered to be fixed since it is assumed that they prefer to trade always with fixed currency of USD or BRICS. A further development of this approach with nodes of fixed opinion is done in [29] for Wikipedia networks when initially there are two groups of nodes of fixed red and blue opinions but all other nodes have neutral state (white opinion or spin zero) and their fixed opinion emerges in the result of asynchronous Monte Carlo process. This process is viewed as the Ising Network Opinion Formation (INOF) model. In fact this INOF model is directly suitable for the modeling of fibrosis progression in the PPI MetaCore network and we are using it here with certain slight modifications described in the next Section. We call this modified system as the Ising Network FIbrosis (INFI) model. It should be pointed out that a somewhat similar asynchronous Monte Carlo process is used in models of associative memory (see e.g., [30,31]).

In this work, using the results reported in [9,16], we take in the INFI model an initial state of 10 fixed red nodes that produce positive activation response of fibrosis on growth factor-$\beta$ (TGF-$\beta$) considering them as fixed 10 red protein nodes. We also consider 6 (or 10, 14) fixed blue protein nodes with negative response (see details in next Section). As a result of fibrosis progression via asynchronous Monte Carlo iterations we obtain the steady-state of affected and healthy protein nodes with about 60–80 percents of proteins affected by fibrosis. This corresponds to a strong fibrosis activation resolution of a large part of the system.

To reduce this strong fibrosis activation we develop an efficient barrage approach based on the analysis of Erdös proteins linked with fixed red nodes. We show that this Erdös barrage constructed in a clever manner allows to obtain a striking reduction of fibrosis activated resolution nodes by a factor of 300 so that almost all protein nodes are transformed in the healthy phase. We hope that the barrage strategy developed in this work will find applications in clinical experiments with fibrosis progression.

The paper constructing as follows: Section 2 describes the data sets, construction of Markov chain transitions, Google matrix and the INFI model, Section 3 presents the obtained results and Discussion and conclusion are given in Section 4. In the Appendix A we provide additional data and Figures.

2. Data Sets and INFI Model Description

2.1. Network Data Sets

In this work we use the same global PPI MetaCore network as in [16]. It contains $N=40,079$ nodes with $N_{\ell }=292,191$ links (without self connections which existed in [15]). The number of directed activation/inhibition links is $N_{\ell +}/N_{\ell -}=65,157/49,321\simeq 1.3$ and the number of neutral directed links is $N_{\ell n}=N-N_{\ell +}-N_{\ell -}=177,713$. Here we do not take into account the bi-functional activation/inhibition nature of links. Thus we simply have a directed network with N nodes and number of directed links being $N_{\ell }$.

For convenience, we present in

Table 1. Table of the subset of $N_r=54$ selected fibrosis proteins (nodes). Here $K_g$ represents the global index of this group, $K_{t,u,d,x}$ represent the index of the four subgroups of 4 TGF-$\beta$ proteins, 20 up-proteins, 20 down-proteins and 10 additional X-proteins; the K ($K^{*}$) index represents the PageRank (CheiRank) index for the global MetaCore network of $N=40,079$ nodes; the last column gives the associated protein names.

${\mathit{K}}_{\mathit{g}}$	${\mathit{K}}_{\mathit{t},\mathit{u},\mathit{d},\mathit{x}}$	K	${\mathit{K}}^{*}$	Protein
1	$K_t=1$	10,780	26,299	TGF-$\beta$ 0
2	$K_t=2$	235	5690	TGF-$\beta$ 1
3	$K_t=3$	968	25,073	TGF-$\beta$ 2
4	$K_t=4$	4726	29,508	TGF-$\beta$ 3
5	$K_u=1$	28,737	25,928	ADAMTS16
6	$K_u=2$	3478	25,137	FGF21
7	$K_u=3$	40,048	28,152	TNFSF18
8	$K_u=4$	2467	19,160	ACAN
9	$K_u=5$	1489	24,511	RPH3A
10	$K_u=6$	26,600	29,559	ADAMTS8
11	$K_u=7$	34,769	39,960	MEGF6
12	$K_u=8$	26,295	27,326	SV2B
13	$K_u=9$	27,111	36,021	C1QTNF3
14	$K_u=10$	34,616	39,841	ANO4
15	$K_u=11$	12,696	16,566	IL11
16	$K_u=12$	26,624	23,640	CDH10
17	$K_u=13$	7263	30,243	HTR2B
18	$K_u=14$	4647	6551	LAMA1-1
19	$K_u=15$	8342	26,295	LAMA1-2
20	$K_u=16$	4021	8252	RAPGEF4
21	$K_u=17$	29,945	36,964	DNER
22	$K_u=18$	22,159	8569	GALNT3
23	$K_u=19$	29,145	15,531	ACSBG1
24	$K_u=20$	24,786	8735	OLFM2
25	$K_d=1$	19,039	28,262	CLEC3B
26	$K_d=2$	26,477	28,290	SCARA5
27	$K_d=3$	26,109	11,185	SLC10A6
28	$K_d=4$	6360	29,204	CXCL5
29	$K_d=5$	14,952	8729	MYOC
30	$K_d=6$	5961	22,288	IFITM1
31	$K_d=7$	5599	4483	ANGPTL4
32	$K_d=8$	25,538	17,434	SELENBP1
33	$K_d=9$	18,938	33,179	FMO1
34	$K_d=10$	34,080	39,427	GPR88
35	$K_d=11$	6276	22,141	HMGCS2
36	$K_d=12$	37,060	28,328	LGI2
37	$K_d=13$	9162	2485	PTN
38	$K_d=14$	513	5974	ADORA2A
39	$K_d=15$	7789	22,652	GFRA1
40	$K_d=16$	6718	8844	IL1R2-1
41	$K_d=17$	35,446	28,306	IL1R2-2
42	$K_d=18$	12,148	3444	PEG10
43	$K_d=19$	27,829	36,195	FMO2
44	$K_d=20$	1973	24,994	COX4I2
45	$K_x=1$	3	13	$\beta$-catenin
46	$K_x=2$	4	6	p53
47	$K_x=3$	11	10	ESR1
48	$K_x=4$	13	25	STAT3
49	$K_x=5$	22	11	RelA
50	$K_x=6$	38	82	PPAR-$\gamma$
51	$K_x=7$	111	767	IKK-$\beta$
52	$K_x=8$	179	198	SNAIL1
53	$K_x=9$	237	1520	MMP-14
54	$K_x=10$	578	2123	Flotillin-1

, taken from [16], 54 selected fibrosis proteins (nodes). These nodes are composed with 4 TGF-$\beta$ proteins/nodes ($K_t=1,2,3,4$), 20 “up-proteins” ($K_u=1,\dots ,20$), 20 “down-proteins” ($K_d=1,\dots ,20$), both obtained from experiments [9] (as described above) and 10 new “X-proteins” (or “X-nodes”; $K_x=1,\dots ,10$) which have a significant influence on the other 44 nodes. The 4 TGF-$\beta$ nodes correspond to different isoforms of this protein. As in [16] we have 4 groups of proteins and we use a specific index for each group: TGF-$\beta$ proteins with index $K_t=1,2,3,4$; up-proteins with a strongest positive response noted by index $K_u=1,\dots ,20$ (ordered by the positive response with the strongest response for $K_u=1$); down-proteins with a strongest negative response noted by index $K_d=1,\dots ,20$ (ordered by the modulus of negative response with the strongest response modulus for $K_d=1$); external proteins noted by index $K_x$ ordered by their local PageRank index (strongest PageRank probability of these 10 proteins is at $K_x=1$; see more detail below). All these 54 proteins have their global index $K_g=1,\dots ,54$ used in Table 1.

2.2. Without Formulas: Methods, Characteristics and Expected Network Results

As in [16] here we present qualitative explanations without formulas of the mathematical methods and characteristics described in the next Subsections. Our aim here is to give a global view on our approach for a common reader.

We use the MetaCore directed network [12] which represents an action of a protein A on a protein B in a form of a directed link (edge) for $N=40,079$ proteins forming the network nodes (proteins). These links are obtained on the basis of careful and detailed analysis of scientific literature about thousands of experiments of various research groups that allowed to collect information about PPI and thus generated a network database with $N=40,079$ nodes and $N_{\ell }=292,191$ links. Certain medical applications of the MetaCore network can be found in [13,14].

The universal mathematical methods to analyze such networks are generic and based on the concept of Markov chains and Google matrix [17,18,19]. The validity of these methods has been confirmed for various directed networks from various fields of science. Therefore, since the Google matrix analysis is based on a generic mathematical foundation, we expect that this analysis will also work efficiently for PPI networks.

The Google matrix of the global MetaCore PPI network G is constructed with specific rules briefly described in the next Section 2.3 and in detail in [17,18,19]. The matrix G is obtained from a matrix of Markov chain transition elements $S_{ij}$ that give weights of transitions between nodes. The important property of G is that its application (multiplication) to an arbitrary initial vector v preserves the probability and the normalization of this vector (sum of all vector elements) remains constant (taken to be unity). As a result of multiple multiplications by G any initial vector converges to the steady-state distribution given by the PageRank vector P. The components of this vector represent the probabilities of each node (protein) in this limit. The nodes with the highest probabilities are the most influential nodes of the network (all nodes are monotonically ordered by decreasing values of the PageRank components which provides the “PageRank index” K such $K\left(j\right)=1,2,\dots$ for nodes j with largest values $P\left(j\right)$). These nodes have typically many ingoing links and it is likely that some of these ingoing links come from other nodes that also have large PageRank values.

For the inverse network, in which all link directions are inverted, the corresponding PageRank is called CheiRank vector $P^{*}$ [19]. The highest probabilities $P^{*}\left(j\right)$ are for nodes j with the CheiRank index $K^{*}\left(j\right)=1,2,\dots$ being the most communicative nodes. They typically have many outgoing links.

In the INFI model Ising spins with polarization $\pm 1$ (red or blue) or 0 (white) are placed at each node of network. They describe fibrosis activated protein state (red or ${\sigma }_j=+1$) which can activate other proteins, healthy protein state (blue or ${\sigma }_j=-1$) which can prepare other proteins and neutral state proteins (white or ${\sigma }_j=0$) that cannot affect other proteins but can be activated or repaired by red or blue proteins. The initial spin configuration has certain fixed red proteins with index $K_g=1,\dots ,10$ (see Table 1) and fixed blue proteins with index $K_g=25,\dots ,30$ (or 34 or 38) whose colors always remain unchanged, all other proteins are white (or ${\sigma }_j=0$). These white proteins, and also red or blue non-fixed proteins, can change their color depending on the majority color of nodes directly linked to it (taking into account the weights of ingoing and outgoing links of the matrix $S_{ij}$ taken without dangling nodes). This choice is similar to social relations when a member of society takes his opinion to be imposed by opinions of other members related with him. An asynchronous Monte Carlo process performs multiple changes of spins leading to the convergence to a steady-state polarization of all spins. The fraction of red spins $f_r$ determines if the fibrosis progression affected a major part of network or not. We give the mathematical definitions of the INFI model in next Subsections.

2.3. Markov Chains, Google Matrix, PageRank and CheiRank

In this section we remind some basic mathematical definitions of Markov chains and how the Google matrix and the related PageRank and CheiRank vectors for the MetaCore network are constructed. For this we use a simple version of the MetaCore network where initial links are not weighted (see Refs. [15,16] for some details about a more refined version with bi-functional links). First, one introduces the adjacency matrix with elements $A_{ij}=1$ if a node j points to another node i (different from j) and $A_{ij}=0$ if there is no link from j to i or if $i=j$ (we use the same convention as in Ref. [16] where all diagonal elements are chosen to be $A_{jj}=0$ even if there is a link from a node j to itself). Then we define a matrix $\tilde{S}$ by ${\tilde{S}}_{ij}=A_{ij}/{\sum }_l{A}_{lj}$ if ${\sum }_l{A}_{lj}\ne 0$ and ${\tilde{S}}_{ij}=0$ otherwise (i.e., for columns j where all elements $A_{lj}$ are 0). This matrix does not yet describe a Markov process since it has potentially some zero columns j (such nodes j are also called dangling nodes). However, in this work we will use a Monte Carlo process based on this particular matrix (see next section).

To obtain a proper stochastic matrix S one replaces $S_{ij}=1/N$ for the dangling nodes j and $S_{ij}={\tilde{S}}_{ij}$ for the other nodes (dangling nodes have no outgoing links so that a column of dangling node j has only zero elements in the matrix ${\tilde{S}}_{ij}$). Here $S_{ij}$ represents the transition probability for a random surfer from node j to i and the column sum normalization ${\sum }_i{S}_{ij}=1$ ensures the conservation of probability.

The Google matrix elements $G_{ij}$ are typically defined by

$$G_{ij}=\alpha S_{ij}+(1-\alpha )/N$$

(1)

where $\alpha =0.85$ is the usual damping factor [17,18]. The Google matrix is also a proper stochastic matrix (column sum normalized) and here the random surfer jumps with probability $\alpha$ on the network according to S and with a probability $(1-\alpha )$ to an arbitrary random node of the network. The damping factor modification helps to avoid possible isolated communities and ensures that the Markov process converges for long times rather quickly to a uniform stationary probability distribution. The latter defines the PageRank vector P which is actually the right eigenvector of the Google matrix G corresponding to the leading eigenvalue $\lambda =1$, i.e., $GP=P$. The elements $P\left(j\right)$ of the PageRank vector correspond to the probability to find the random surfer on the node j in the stationary limit of the Markov process. The PageRank index $K\left(j\right)$ is obtained by ordering the nodes with decreasing values of $P\left(j\right)$, i.e., the highest (lowest) PageRank probability $P\left(j\right)$ corresponds to $K\left(j\right)=1$ ($K\left(j\right)=N$). In this work we will use the K-rank of a node to characterize this node in a unique way. The PageRank probability $P\left(j\right)$ is typically related to the number of ingoing links pointing to node j. However, it also takes into account the “importance” (i.e., PageRank probability) of the nodes having a direct link to j.

One can also consider the network obtained by the inversion of all link directions. The same construction as described above provides a Google matrix noted as $G^{*}$ and the corresponding PageRank vector is called the CheiRank vector $P^{*}$, defined by $G^{*}{P}^{*}=P^{*}$. The CheiRank probability $P^{*}\left(j\right)$ is typically related to the number of outgoing links weighted by the value $P^{*}\left(i\right)$ of nodes i having a link from j to i. The CheiRank index $K^{*}\left(j\right)$ is obtained by ordering the CheiRank vector with decreasing values of $P^{*}\left(j\right)$.

For further details about the properties of Google matrix, PageRank and CheiRank vectors, with many different example networks (also with proteins networks), we refer to [15,16,18,19].

2.4. Ising Spin Network, Monte Carlo Process for INFI Model

In this work, we consider an Ising type of model using the symmetrized matrix $\tilde{S}+{\tilde{S}}^T$ to model the effective spin interactions. We denote by ${\sigma }_i$ the spin associated to the node (protein) i but here we allow for three different values being ${\sigma }_i=+1$ (“red state”, “fibrosis state”), ${\sigma }_i=-1$ (“blue state”, “healthy state”) and also $\sigma =0$ (“white state”, “undetermined or neutral state”). We study this model using a specific Monte Carlo process similar to [29]. The white state will only be used in the initial condition and is supposed to be only temporary.

First, we fix 10 specific nodes permanently to the red state, those with $K_g=1,\dots ,K_g=10$ in the set of Table 1 corresponding to the four TGF proteins and the first 6 $K_u$-proteins with strongest positive response on the TGF stimulation [16]. We also fix $n_b$ nodes permanently to the blue state being the nodes with $K_g=25,\dots ,K_g=24+n_b$ corresponding to the first $n_b$ $K_d$-proteins with strongest negative TGF-response. Here we choose mostly $n_b=6$ (same number as the 6 $K_u$ proteins with fixed red state) but we also present some data for the cases $n_g=10$ (same total number of fixed red proteins) and $n_g=14$. Therefore $10+n_b$ nodes have already a permanently fixed spin value ${\sigma }_i=\pm 1$.

We implement a Monte Carlo procedure that modifies the other $N_v=N-10-n_b$ remaining nodes with “variable” spin values (non fixed). For this at given iteration time $\tau$, we update the $N_v$ variable nodes by computing the sum over all nodes j of the network (fixed and variable):

$$Z_i=\sum _{j(\ne i)}({\tilde{S}}_{ij}+{\tilde{S}}_{ji}){\sigma }_j$$

(2)

and choosing ${\sigma }_i(\tau +1)=+1$ if $Z_i>0$, ${\sigma }_i(\tau +1)=-1$ if $Z_i<0$ and ${\sigma }_i(\tau +1)={\sigma }_i\left(\tau \right)$ (unchanged) if $Z_i=0$. Without going into technical details, we mention that the numerical implementation of (2) exploits of course the sparse structure of the matrix ${\tilde{S}}_{ij}$. This update is done for every node i in the set of $N_v$ variable nodes in a random order using a random permutation (thus no repetitions) which is chosen at the beginning of the iteration process. Once a spin value for a node i has been updated, the new value ${\sigma }_i(\tau +1)$ will be used for the update of the subsequent nodes $\tilde{i}$. We note that somewhat similar asynchronous Monte Carlo process what used for bio-networks in [32,33] but there the matrix elements ${\tilde{S}}_{ij}+{\tilde{S}}_{ji}$ had values $\pm 1$ and network sizes were relatively small (about 100).

This full update procedure for one time step and the full set of $N_v$ variable nodes is repeated up to the iteration time ${\tau }_{\max }=100$ (with different random permutations for the update order at each time step). We also verify if at any given value of $\tau$ there is ideal convergence where the spins no longer change. (This happens typically at ${\tau }_{last}\approx 10$; see below for details).

The full iteration procedure with ${\tau }_{\max }=100$ time steps is repeated $R=100,000$ times with different random permutations and eventually also different random positions of initial blue or white nodes (depending on the precise type of initial condition).

During this procedure, we compute the fraction of red outcome $f_r\left(i\right)$ for each node as the number of times the node i has red state at ${\tau }_{\max }$ (at different realizations) divided over R and also the overall average $f_r={\sum }_i{f}_r\left(i\right)/N_v$ over the nearly full network of variable nodes. To test the convergence, we also compute the overall average at intermediate values of $\tau$.

In order to test the procedure and also in the direct influence of the permanently fixed nodes, we first choose an initial condition where all $N_v$ variable nodes are initialized to the white state. Figure A1 illustrates for this case the convergence of (the overall network average) $f_r\left(\tau \right)$ with iteration time $\tau$. The convergence seems to be very good at $\tau >5$ with nearly constant values on graphical precision. However, a closer look at the data shows that $|f_r(\tau +1)-f_r\left(\tau \right)|\sim {10}^{-5}$ for $\tau =10$ and the value of $f_r\left(\tau \right)$ becomes constant only at $\tau \ge {\tau }_{last}\approx 70$ with typical last non-vanishing differences $|f_r({\tau }_{last}-1)-f\left({\tau }_{last}\right)|\sim {10}^{-10}$ and $|f_r({\tau }_{last}+1)-f\left({\tau }_{last}\right)|=0$. We have verified that this behavior is due to the fact that typically at $\tau \approx 10$ there is exact convergence for a given realization of the random permutations (random pathways) with no further modifications of the spins at $\tau >10$. However, there are rare realizations with larger values, up to $\tau ={\tau }_{last}\approx 70$ for exact convergence.

We stress that the approach of INFI model described above is based on physics of interacting Ising spins placed at each node of MetaCore network so that the configuration space of the system becomes $2^N$ instead of size N used in the Google matrix methods. Thus the INFI model provides a qualitatively new approach for analysis of PPI networks. For Wikipedia networks a similar approach was used in [29]. We state that this method previously was never used for analysis of PPI networks.

In the next section, we present the results for different quantities and also other initial conditions.

3. Results

3.1. Numerical Results

All results presented in this work were obtained with maximal iteration time ${\tau }_{max}=100$ and $R=100,000$ realizations for (potential) different initial conditions and random permutations for the Monte Carlo procedure. Also when we describe a particular initial condition it is implicitly understood that the permanently fixed 10 red and $n_b$ blue nodes given above are indeed fixed to these values from the very beginning, i.e., if we say “that there are $n_{ib}$ initial blue values” (with $n_{ib}=0,1,\dots$) we mean that there are $n_{ib}$ blue nodes on the set of $N_v$ variable (non-fixed) nodes, other variable nodes are initialized to white values and the fixed nodes have still their 10 red and $n_b$ blue values (in particular the value of $n_{ib}$ does not include the number of $n_b$ fixed blue nodes).

In Figure 1, we show for the case of white or zero initial spins (for the $N_v$ variable nodes; i.e., $n_{ib}=0$) and the three cases of $n_b=6,10,14$ initial blue fixed nodes (and 10 fixed nodes as explained above) the probability distribution (normalized by an integral) of fraction of red outcomes $f_r\left(i\right)$ for the nodes obtained by a histogram with bin width $\Delta f_r=0.01$. These distributions are strongly peaked at values close to $f_{r,peak}\approx 0.93$ ($n_b=6$), $f_{r,peak}\approx 0.81$ ($n_b=10$), $f_{r,peak}\approx 0.65$ ($n_b=14$) for roughly 82% of nodes. The secondary peak at $f_r=0$ corresponds to the fraction (number) $\approx 0.18$ (7230) of nodes (same value for the three cases of $n_b$) which stay white after 100 iterations for all $R=100,000$ pathway realizations. This set of stable white nodes correspond to nodes not connected to the small number of fixed red and blue nodes. However, despite the small value of only 10 initial fixed red notes the majority of the other nodes has a red outcome, especially for $n_b=6$ and with slightly reduced values for $n_b=10,14$. Note that the bin width $\Delta f_r=0.01$ in Figure 1 is still the rather large and histogram computations with $\Delta f_r=0.001$ and $\Delta f_r=0.0001$ show that the actual peaks of the distributions are much sharper as visible in Figure 1.

It seems that with the white initial condition (except for the fixed nodes) nearly all nodes (except those in the stable white set) have a large probability for red outcome. Therefore, we try to reduce this red outcome by choosing a certain number $n_{ib}$ of initial blue (variable) nodes and other nodes with white initial values. The question is also where to place these blue nodes. In a first model we choose the $n_{ib}$ blue nodes at different random positions (for each of the $R=100,000$ pathway realizations) on the full set of variable nodes giving a set of $N_r=N_v$ of possible blue initial nodes with possible values $n_{ib}=0,\dots ,N_r$ and a corresponding fraction $f_{ib}=n_{ib}/N_r$. Figure 2 shows for the case $n_b=6$ (the results for $n_b=10$ and $n_b=14$ are very similar) the overall network average $f_r$ at maximum iteration time as a function of $f_{ib}$ (or $n_{ib}$). The initial value $f_r\approx 0.75855$ at $n_{ib}=0$ corresponds to the average which can also be obtained from the histogram data shown in Figure 1 which is roughly $0\times 0.18+0.93\times 0.82\approx 0.76$. With increasing $f_{ib}$ (or $n_{ib}$) the value of $f_r$ decreases, e.g., $f_r\approx 0.1$ for $n_{ib}\approx 60$ and $f_r\approx {10}^{-2}$ for $n_{ib}\approx 100$. However, this decrease is not optimal since the initial $n_{ib}$ blue nodes can be on arbitrary random positions of available variable nodes. Therefore, we also try different network subsets for the potential initial blue nodes.

In Figure 3, we show the dependence of $f_r$ on $f_{ib}$ (or $n_{ib}$) for a reduced subset of $N_r=38$ nodes corresponding to the 38 proteins of Table 1 which are not used for fixed red/blue values, i.e., $K_g=11,\dots ,K_g=24$ and $K_g=31,\dots ,K_g=54$. Now, $n_{ib}$ represents the number of initial blue nodes on random positions in this subset (and white initial nodes on every other variable node). Now, the decrease of $f_r$ with $n_{ib}$ seems somewhat stronger, i.e., $f_r\approx 0.1$ for $n_{ib}\approx 8$ and $f_r\approx {10}^{-2}$ for $n_{ib}\approx 16$, and we may conclude that this subset is more effective to reduce the red outcome (“to block fibrosis”) than the full set of available variable nodes.

To determine a still more effective subset, we determine all nodes which have a direct link or inverse link to one of the 10 fixed red nodes (those with $K_g=1,\dots ,K_g=10$ in Table 1). This defines a particular set, of size $N_E=353$, with Erdös number being unity with respect to the 10 fixed red nodes as HUB and using the symmetrized link matrix ${\tilde{S}}_{ij}+{\tilde{S}}_{ji}$. Figure 4 shows the dependence of $f_r$ on $f_{ib}$ (or $n_{ib}$) for this Erdös set with $N_r=N_E$. Now the decrease is even more effective with $f_r\approx 0.1$ for $n_{ib}\approx 6$ and $f_r\approx {10}^{-2}$ for $n_{ib}\approx 14$.

In Figure 5, we show in a color plot the dependence of $f_r\left(i\right)$ for the 54 nodes i belonging to the set of Table 1 on the index $n_g$ which is a monotonic function of $n_{ib}$ (essentially linear for $n_{ib}\le 10$ and logarithmic for $n_{ib}>10$; see caption of Figure 5 and Appendix Figure A1 for details). The data of Figure 5 correspond to the data of Figure 4 using the Erdös set with $N_E=353$ as potential initial blue nodes (i.e., for each value of $n_{ib}\in \{0,\dots ,353\}$ we have $n_{ib}$ initial blue nodes at random positions in this set). In Figure 5, we can of course identify the 10 fixed red nodes ($K_g=1,\dots ,10$) and the $n_b=6$ fixed blue nodes ($K_g=25,\dots ,30$) with either $f_r\left(i\right)=1$ or $f_r\left(i\right)=0$ respectively for all values of $n_g$. The other nodes follow quite closely the decrease of Figure 4 for the global average of $f_r$. However, for certain specific nodes ($K_g=17,23,36,38,41,44$) the decay of $f_r$ with $n_g$ (or $n_{ib}$) is less pronounced for $n_g>10$. Apparently, these nodes are less likely to be blocked by a modest number of initial blue nodes.

We have also analyzed the data of Figure 2 (using all $N_v$ variable nodes as potential initial blue nodes) and Figure 3 (using the remaining set of 38 non-fixed nodes of Table 1 as potential initial blue nodes) with similar color plots and in both cases we observe the same qualitative behavior as in Figure 5: identification of fixed red/blue nodes, similar decrease of $f_r\left(i\right)$ with increasing $n_{ib}$ for the other nodes and less pronounced decrease for the 6 specific nodes mentioned above.

The question arises which of the nodes of the Erdös set, or more generally, which configurations of few selected nodes of this set, are most effective to reduce the red outcome when selected as initial blue nodes. To answer this question, we compute for each node i of the Erdös set the red outcome $f_{rc}\left(i\right)$ averaged over the full network when this specific node i is selected as single initial blue node, i.e., with $n_{ib}=1$ but now with different cases of given fixed positions (instead of random positions). Note that this quantity is different from $f_r\left(i\right)$ used in the histogram of Figure 1 which is the probability of red outcome of node i (not averaged over the network) with full white initial condition. The nodes of the Erdös set can be ordered with increasing values of the new quantity $f_{rc}\left(i\right)$ which provides a specific ranking index $K_{f_r}$ in this set. In

Table 2. Table of selected fibrosis proteins (nodes) belonging to a smaller subset of the Erdös group such that either $K\le 40$ or $K^{*}\le 40$ where the K ($K^{*}$) index represents the PageRank (CheiRank) index for the global MetaCore network of $N=40,079$ nodes. $f_{rc}$ represents the average fraction of red nodes obtained when using the corresponding node as single initial non-fixed blue node and the index $K_{f_r}$ is obtained by ordering the Erdös group of 353 elements with increasing values of $f_{rc}$; the last column gives the associated protein names. The shown subset of the Erdös group corresponds essentially to the nodes with lowest $f_{rc}$-value (up to $K_{f_r}=17$). The four nodes marked by “*” in the first column are the nodes used for the subsequent example computations using 4 specific optimal nodes.

	${\mathit{K}}_{{\mathit{f}}_{\mathit{r}}}$	K	${\mathit{K}}^{*}$	${\mathit{f}}_{\mathbf{rc}}$	Protein
*	1	17	1	0.010691	c-Myc
*	2	4	6	0.029542	p53
*	3	342	7	0.038899	c-Fos
	4	14	14	0.053526	Androgen receptor
	5	22	11	0.075088	RelA (p65 NF-kB subunit)
	6	94	34	0.07674	HDAC1
	7	188	19	0.078157	p300
	8	6272	3	0.078989	IGF2BP3
	9	232	9	0.081749	SP1
	10	64	29	0.083725	HIF1A
*	11	3	13	0.083939	$\beta$-catenin
	12	203	21	0.13861	E2F1
	13	728	5	0.14694	SOX9
	14	432	24	0.14808	BRG1
	15	8	115	0.15027	EGFR
	16	72	18	0.15149	EZH2
	17	13	25	0.15649	STAT3
	20	480	40	0.17501	C/EBP $\beta$
	21	38	82	0.20211	PPAR-$\gamma$
	24	298	33	0.20997	ELAVL1 (HuR)
	26	394	37	0.22213	CREB1
	31	6370	17	0.28502	PUM2
	32	17,711	8	0.29183	CUX1 (p110)

, we present in $K_{f_r}$-order the nodes (proteins) of the Erdös set which have either $K\le 40$ or $K^{*}\le 40$, either low K- or $K^{*}$-rank. It turns out these nodes also provide the lowest values of $f_{rc}$ and $K_{f_r}$ (with some holes for $K_{f_r}>17$).

Then we choose as example four “optimal” nodes with $K=3,4,17,342$ (marked with an asterix in the first column of Table 2). The nodes $K=4,17,342$ occupy indeed the top three places in the $K_{f_r}$-rank while the node $K=3$ “only” corresponds to the position $K_{f_r}=11$. The reason for this choice is related to the fact that these four nodes are more uniformly optimal if we also choose a small number of $n_{ib}=2,\dots ,10$ initial blue nodes. In this case, using the data of Figure 4 (or the code to produce these data) it is actually possible to compute the conditional probability $f_{rc}\left(i\right)$ of a red outcome when the node i is by chance selected by the random initial condition for $n_{ib}>1$ (as one of the initial $n_{ib}$ blue nodes). For example for $n_{ib}=4$ with $R=100,000$ different random initial conditions of 4 blues nodes out of 353 we have typically a bit more than 1000 realizations where an arbitrary fixed node i belongs to the random set of 4 initial nodes. This provides enough data for a reasonable average to compute the conditional probability of red outcome of the node i. In this way, it is possible to compute more general $K_{f_r}$-rankings as in Table 2, also for modest values of $n_{ib}>1$. It turns out that these rankings produce roughly the same sets of nodes in the first places (with possible permutations between different $n_{ib}$-values) and the four selected nodes $K=3,4,17,342$ are indeed optimal as a group for the three values $n_{ib}=2,3,4$ (and also some larger values).

In Figure 6, we present (for $n_b=6$) the dependence of $f_r$ on $n_i$ for this small optimal set of 4 possible initial blue nodes (red square data points). Here the decrease of $f_r$ with increasing $n_i$ is indeed very strong with $f_r\approx 0.035$ for $n_{ib}=1$ and $f_r\approx 0.00318$ for $n_{ib}=4$. Furthermore, Figure 6 also provides the conditional probabilities of red outcome $f_{rc}([K_1,\dots ,K_{n_{ib}}]$ (black small circle data points) for specific configurations $[K_1,\dots ,K_{n_{ib}}]$ with one configuration for $n_{ib}=0,4$, four configurations for $n_{ib}=1,3$ and six configurations for $n_{ib}=2$. We see that there are considerable fluctuations between the different configurations in the red outcome for $n_{ib}=1$ and $n_{ib}=2$ with optimal values being $f_{rc}\left(\left[17\right]\right)\approx 0.009477$ and $f_{rc}\left([3,17]\right)\approx 0.004863$. Appendix Figure A3 provides two panels of similar figures for the other cases $n_b=10$ and $n_b=14$.

Figure 7 is a similar color plot as Figure 5 but for the specific set of 4 optimal nodes $[3,4,17,342]$. The conclusions are similar to Figure 5 (confirmation of fixed red/blue nodes, strong decrease with increasing $n_{ib}$ for other nodes and certain specific nodes with less pronounced decrease.

In Figure 8, we show the 353 nodes of the Erdös set in the global K-$K^{*}$-plane (in a double logarithmic representation) with colored data points such that the color provides the value of $f_{rc}\left(K\right)$ for each node at given K-value (obtained at $n_{ib}=1$ with the corresponding node as initial blue node). Optimal nodes with low red outcome (blue color) have typically small values in their K- and $K^{*}$-rank and nodes with large red outcome (red color) have typically large values in their K- and $K^{*}$-rank. Furthermore, Figure 8 also provides the positions of the 10 fixed red and the 6 blue nodes in the global K-$K^{*}$-plane with typically quite large values for K and $K^{*}$. Using the same data Appendix Figure A4 shows the dependence of $f_{rc}$ (at $n_{ib}=1$) on the rank $K_{f_r}$. Obviously, this curve is monotonically increasing and the range with small values of $f_{rc}$ is rather small, i.e., $f_{rc}<0.1$ corresponds to $K_{f_r}\le 11$ (see also Table 2). Note that the average of the curve in Figure A4 corresponds to the data point at $n_{ib}=1$ in Figure 4 which is $f_r(n_{ib}=1)\approx 0.5818$.

We have also analyzed the statistical properties of the individual $f_r\left(i\right)$ values, same quantity as in Figure 4, but using the data of Figure 6 for $n_{ib}=4$ with the initial blue node configuration $[3,4,17,342]$ and global (network averaged) value $f_r\approx 0.00318$. For example Appendix Figure A5 shows the 100 nodes with largest values of this quantity in the K-$K^{*}$ plane. Obviously, among the 26 nodes with $f_r\left(i\right)=1$ we have the 10 permanently fixed red nodes ($K_g=1,\dots ,10$ in Table 1) but there are also further 16 nodes with $f_r\left(i\right)=1$ probably with (mostly) exclusive links to the fixed 10 red nodes, i.e., no or few links to other nodes, therefore explaining the fixed outcome $f_r\left(i\right)=1$. There are about 1% (4%) of network nodes with values $f_r\left(i\right)>0.5$ ($f_r\left(i\right)>0.08850$).

Furthermore, Appendix Figure A6 shows the histogram distribution $p\left(f_r\right)$ for the same data with a rapid decay at $f_i<1$ and a small peak at $f_r=1$ corresponding to the 26 nodes i with $f_r\left(i\right)=1$. We have also computed the related global probability $P\left(f_r\right)$ for a node i to have a value $f_r\left(i\right)>f_r$. This quantity is shown in Appendix Figure A7, confirming the rapid decay which is quite well algebraic as $P\left(f_r\right)\sim f_r^{-1.5}$ (for $f_r\ge 5\times {10}^{-3}$) and which is similar to the Poincare recurrences decay in symplectic chaotic maps with Ulam networks [34].

The slow algebraic decay of $P\left(f_r\right)$ has important consequences. It shows that the average values of $f_r$ shown in Figure 6 even underestimate the effect. Explicitely, the average value of $f_r\approx 0.00318$ for $n_{ib}=4$ in Figure 6 corresponds, according to Appendix Figure A7, to a fraction of nodes $P\left(0.00318\right)\approx 0.22$ having an $f_r\left(i\right)$ value above this average and 78% of nodes have a smaller $f_r\left(i\right)$ value. It is well known, that for such long tail distributions one should also focus on the median value $f_{r,median}$ defined by $P\left(f_{r,median}\right)=0.5$ corresponding in our case to $f_{r,median}={10}^{-5}$ (only 1 red outcome in the $R=100,000$ pathway realizations) indicating that (slightly more than) 50% of nodes have an $f_r\left(i\right)$ value below or equal to the median value with nearly perfect reduction of the red outcome.

3.2. Results Without Formulas

In this work we presented a mathematical model of fibrosis progression in the PPI MetaCore network describing the global interaction structure of almost all proteins and important molecules (nodes). We show that even with only 10 fibrosis activated proteins the fibrosis progression can spread over a great majority of nodes (about 70%). The developed analysis of network structure allows us to propose an efficient strategy which allows to reduce a number of fibrosis activated nodes by a factor 300 and disease elimination. The method is based on the Erdös barrage construction: we determine the Erdös nodes directly linked to the fixed 10 activated red nodes; this number can be relatively large (353 in our case); however, we show that a barrage with only 4 blue repairing Erdös nodes, corresponding to the four proteins c-Myc, p53, c-Fos and $\beta$-catenin, gives a reduction of the average number of fibrosis activated nodes by a factor 300; these 4 nodes belong to network nodes with high PageRank and CheiRank indexes of global MetaCore network.

Furthermore, this average actually underestimates the effect since it is determined by a relative small number of nodes with a modest reduction (e.g., factor of $\sim 100$) while for more than 50% of nodes the reduction factor is even 100,000 (only 1 infection outcome in the 100,000 statistical realizations of our simulation). We have also identified in the group of Table 1 six interesting proteins HTR2B, ACSBG1, LGI2, ADORA2A, IL1R2-2 and COX4I2 (corresponding to the vertical green lines in Figure 5 and Figure 7) where the reduction effect is somewhat less pronounced (compared to the 50% of nodes with nearly perfect reduction). We expect that our INFI model can be tested with other PPI networks (e.g., those of [35,36]).

4. Discussion and Conclusions

Decoding fibrosis-associated proteins and fibrotic remodeling progression is a critical issue in treating heart failure. An experimental determination of the fibrosis resolution is extremely time consuming and dynamic process. In this work we present and describe a mathematical model of the fibrosis progression using the global PPI MetaCore network. The performed analysis shows that in TGFβ-mediated fibroblast activation, PPI hubs predict new regulatory proteins for fibrosis progression. We developed an efficient method of the Erdos barrage to identify the proteins driving fibrosis activation cascades in a physiopathological scenario using MetaCore network.

In cardiac tissue, resident cardiac fibroblasts respond to stress or injury by evoluting through phases determined by a inflammatory motifs, a proliferative capacity, and finally, a reparative status, characterized by the deposition of ECM that culminates in fibrosis. The molecular mechanisms orchestrating the dynamic control of the fibroblast phenotype remain obscure, hampering the development of therapeutic strategies to combat cardiac fibrosis. In our study, we found that c-Myc, p53, c-Fos and Androgen receptor arbitrate cardiac fibroblast activation in cardiac remodeling. The gene encoding the c-MYC is the most frequently perturbed transcription factor in the majority of human cancers, with its activation of downstream target genes in the vast majority of the cancer types [37]. Dysfunction of c-MYC can occur through elevated transcription of the gene, amplification or protein stabilization, all of which release the otherwise tight control and altered expression levels of the protein in the cell [38]. nhibition of c-MYC has been observed to favor apoptosis, growth arrest, differentiation, senescence, metabolic changes and tumour regression in several human cancer models, clearly reflecting its potential as a target for anti-cancer therapy [39]. In activated cardiac fibroblasts, p53 transcriptional targets are among the most deregulated gene programs in pressure overload-induced tissue remodeling [40]. In cancerous tumors, p53 is an important component of the DNA damage and transcriptional regulator of cell cycle [41]. Interestingly, p53-sensitive cellular senescence has also been linked to myocardial remodeling in animal models ou tissue fibrosis [42]. Indeed, over expression of the p53 target gene in transgenic mice inhibits fibroblast activation and reduces scar formation after myocardial infarction [43]. An important factor C-fos is a member of a family of early genes which are involved in the signal transduction cascade coupling extracellular stimuli to intracellular events [44]. In particular, C-fos is critical to the regulation of transcription, and may mediate long-term effects of growth factors and membrane-depolarising signals on neural responses to stress [45]. Androgen receptor (AR) signalling in fibroblasts is important in carcinogenesis and prostate development [46]. However, the molecular mechanisms of AR action in fibroblasts and other non-cardiacl cell types are largely unknown.

In the past few years, the receptors and signal transduction pathways mediating the effects of TGF-$\beta$ on fibroblasts have been identified, enabling the dissection of the specific pathways involved in pathogenic events. TGF-$\beta$ type I and type II transmembrane receptor serine/threonine kinases transduce downstream signals via novel cytoplasmic latent transcription factors, Smad proteins [47]. Both Smad2 and Smad3 are phosphorylated directly by the type I receptor kinase and able to translocate to the nucleus, where they act as transcriptional regulators of target genes controlling cell apoptosis and differentiation [48]. Interestingly, deletion of Smad3 results primarily in impaired status of immunity in mice and shortening their lifespan [49]. The Smad pathway is believed to be the major signaling mechanism through which active TGF-$\beta$ stimulates the induction of profibrotic cascades.

In a profibrotic scanario, TGF-$\beta$-activated kinase 1 (TAK1) is one of the best characterized non-Smad signal transducers critical for TGF-$\beta$ functions in EMT and apoptosis through activating the c-Jun N-terminal kinase (JNK) and p38 MAPK cascade [50]. TAK1 also plays an essential role in mediating TGF-$\beta$ activation of I-kappa B kinase (IKK) and the master transcription factor nuclear factor kappa B (NF-kB) that is required for cell survival [51]. In analogy to the mechanism defined in interleukin-1/Toll-like receptor pathways, TGF-$\beta$-induced activation of TAK1 requires TRAF6, a RING domain ubiquitin ligase that itself is modified by a K63-linked polyubiquitin chain, which acts as a scaffold to recruit TAK1 to the TGF-$\beta$ receptor complex and triggers TAK1 activation [50]. Activity of TAK1 is also regulated by its binding proteins, including TAK1-binding protein 1 (TAB1) that binds constitutively the kinase domain and TAB2 or TAB3 that binds the C-terminal domain and functions as an adaptor linking TRAF6 to TAK1 [52]. However, how TGF-$\beta$ coordinates TAK1 kinase to dictate the opposing responses of cell survival and apoptosis in different cellular scenarios.

Thus the PPI network leads to a very complex system properties. In view of that we hope that the proposed INFI model may lead to better understading of PPI systems.

The presented method of INFI model is based on mathematical properties of PPI networks and thus it is generic and universal. It can be applied to any other disease progression on PPI network and construction of Erdös barrage for disease propagation. As initial income for the application of this method a medical researcher should provide positive and negative responses on certain protein stimulation similar to those of Table 1 with TGF-$\beta$ stimulation of fibrosis.

Thus we developed an efficient method of the Erdös barrage when a small group of e.g., 4 repairing protein can reduce the number of fibrosis activated proteins (nodes) by a factor 300 leading to a healthy state of the global system described here by the MetaCore network. We expect that similar results can be obtained for other disease progression. We think that it would be very interesting to test this INFI approach with other global PPI networks like [35,36]. We hope that the described INFI method will lead to new efficient medical treatments of fibrosis and other various diseases.