- Email: [email protected]
Statistical behavior of time dynamics evolution of HIV infection
JID:PLA
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.1 (1-5)
Physics Letters A ••• (••••) •••–•••
Contents lists available at ScienceDirect
Physics Letters A www.elsevier.com/locate/pla
Statistical behavior of time dynamics evolution of HIV infection Ramón E.R. González, Iury A.X. Santos, Marcos G.P. Nunes, Viviane M. de Oliveira ∗ , Anderson L.R. Barbosa ∗ Departamento de Física, Universidade Federal Rural de Pernambuco, 52171-900, Recife, Pernambuco, Brazil
a r t i c l e
i n f o
Article history: Received 15 May 2017 Received in revised form 14 July 2017 Accepted 17 July 2017 Available online xxxx Communicated by C.R. Doering Keywords: Random matrix theory Cellular automata HIV infection
a b s t r a c t We use the tools of the random matrix theory (RMT) to investigate the statistical behavior of the evolution of human immunodeficiency virus (HIV) infection. By means of the nearest-neighbor spacing distribution we have identified four distinct regimes of the evolution of HIV infection. We verified that at the beginning of the so-called clinical latency phase the concentration of infected cells grows slowly and evolves in a correlated way. This regime is followed by another one in which the correlation is lost and that in turn leads the system to a regime in which the increase of infected cells is faster and correlated. In the final phase, the one in which acquired immunodeficiency syndrome (AIDS) is stablished, the system presents maximum correlation as demonstrated by GOE distribution. © 2017 Elsevier B.V. All rights reserved.
1. Introduction Since the first cases of human immunodeficiency virus (HIV) infection were diagnosed many investigations have been carried out in order to understand the evolution of the infection [1–6]. The dynamics of HIV infection is characterized by three distinct phases which are related to the evolution of the degree of infectability. In the first one, the primary infection, a peak in the infection level is observed around 5 weeks after the contagion. This phase lasts from 10 to 12 weeks. After this period the immune system presents false signs of recovery. The second phase, the clinical latency period, is characterized by a very low viral concentration followed by a gradual growth which leads to the impairment of essential functions of the immune system. This period may last for 2 to 10 years. After this period the probability of the patient developing acquired immunodeficiency syndrome (AIDS) is very large. Several mathematical models have been proposed to study the dynamics of HIV infection [4]. Most of the investigations are based on the use of differential equations but models in which the spatial structure is considered have also been developed. In this sense, Zorzenon dos Santos and Coutinho [1] used a cellular automata model to study the evolution of HIV infection and the onset of AIDS. Their model propose a local interaction among the different types of cells (healthy, infected and dead), while the lymph node is modeled by using a square lattice with periodic boundary condi-
*
Corresponding authors. E-mail addresses: [email protected] (V.M. de Oliveira), [email protected] (A.L.R. Barbosa). http://dx.doi.org/10.1016/j.physleta.2017.07.022 0375-9601/© 2017 Elsevier B.V. All rights reserved.
tions and Moore neighborhood. The model was able to reproduce the time scales of the three phases of the infection. The random matrix theory (RMT) was proposed by Wigner and Dyson for the study of the statistical properties of the spectrum of complex nuclei [7,8]. This approach has been successfully applied to the study of many subjects such as wireless communications [9], chaotic systems [10], protein dynamics [11,12] and structure and interactions of biological networks [13]. Furthermore, the RMT can also be applied to multivariate time series as electroencephalographic recordings with epileptic events [14] and financial time series [15,16], which motivated us using this approach to the time dynamics of HIV infection. In this work we use the cellular automata model introduced by Zorzenon dos Santos and Coutinho [1] to study the statistic behavior of time dynamics of HIV infection by means of the RMT. In Section 2, we describe the cellular automata model used to study the evolution of HIV infection and we also present the time evolution of density of infected cells as generated by the model. A brief review about equal-time cross-correlation matrix and RMT is made in Section 3. Finally, in Section 4 we show our results and discussion. 2. Cellular automata model for the dynamics of HIV infection In the cellular automata model used by Zorzenon dos Santos and Coutinho [1] the square lattice is composed by cells that can be in one of four possible states: healthy, infected-1, infected-2 and dead. In the infected-1 state, the infected cells are able to spread the infection. The state infected-2 corresponds to infected cells which have already been recognized by the immune system.
JID:PLA
2
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.2 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
In the initial configuration, the system is composed by healthy cells and also by a small fraction of infected cells, p H I V , which are randomly distributed on the lattice. The infection process is governed by four rules: (1) A healthy cell will become an infected-1 cell if there is at least one infected-1 cell in its neighborhood (Moore neighborhood), or if there are at least R infected-2 cells in its neighborhood; (2) After spend τ time steps in the same state, an infected-1 cell will become an infected-2 cell, in which τ is the time necessary for the immune response; (3) An infected-2 cell will become a dead cell in the next time step; (4) A dead cell will become a healthy cell with probability p repl . If it happens, the healthy cell will become an infected-1 cell with probability p inf ec . After the update of the lattice according to the rules described above one time step is counted which corresponds to one week. In Fig. 1 we show the density of infected cells as a function of time, obtained from the model of Zorzenon dos Santos and Coutinho [1], in which the initial conditions are p H I V = 0.05, p inf ec = 10−5 and p repl = 0.99, τ = 4 and R = 4. We have used a lattice of size equal to L = 700 and the results are obtained from averages over 500 distinct simulations. We can clearly notice a qualitative behavior very similar to the one presented by the dynamics of the HIV infection. In the period between 0 and 12 weeks we verify the first phase of the dynamics (primary infection), that is characterized by a peak of the infected cells. In the second phase, the clinical latency, we observe three regions which are delimited in Fig. 1: In region (a) the concentration of infected
cells is very small and presents a very slow increment with time; in the second regime (b) we verify an augment of the fluctuation of the concentration of infected cells in which the curve presents a transition region with an inflection point leading to another regime (c) that is characterized by a faster increasing of the concentration of infected cells leading to the development of AIDS. In the third phase (d) we observe the onset of AIDS. We can better understand the behavior of the system by looking at the snapshots of five distinct time steps of some realizations of the dynamics. In Fig. 2(a) we can visualize the system during the primary infection period. After the introduction of the infected cells, pulses of these cells with width (τ + 1) are generated and propagate in all directions. Around the fifth week of the infection process one can notice a peak in the concentration of infected cells and almost all the sites of the lattice are composed by cells of type infected-1 and infected-2 (which are represented in yellow and green, respectively). From this time the concentration of infected cells decreases and the end of the primary infection phase is observed when the number of time steps is equal to 2(τ + 1). After this initial phase new infected cells will be originated when dead cells are replaced by infected-1 cells. They generate structures like the ones showed in Fig. 2(b), which spread over the lattice in a much longer time that the one observed in primary infection. As long as these structures increase they can occupy the whole lattice (see Figs. 2(c) and 2(d)), incrementing the number of infected cells until they reach a stationary state. At the end of this long period in which the increase of these structures is observed, we verify another pattern in which the distinct cells are randomly distributed in the lattice with the predominance of the infected cells (Fig. 2(e)), characterizing the onset of AIDS. 3. Cross-correlation and random matrix theory In this section, we will do a brief review about equal-time cross-correlation matrix [14] and RMT [8]. Their concepts will be applied to study the time dynamics evolution of HIV infection statistic in the next section. The equal-time cross-correlation matrix C is obtained from an M-dimensional temporal series set H i j with length T , where i = 1, . . . , M and j = 1, . . . , T with M ≤ T . The correlation matrix is given by
C=
1 T
H Ht ,
(1)
where index t is a transposed of matrix and the H matrix elements are built from temporal series set as following
Fig. 1. Time evolution of the density of infected cells obtained from the simulations. The results are obtained from averages over 500 simulations. We used lattices with size L = 700, p H I V = 0.05, p inf ec = 10−5 , p repl = 0.99, τ = 4 and R = 4. PI means primary infection.
Hij =
Hij − Hi
σHi
.
(2)
The H i and σ H i are the average and standard deviation of temporal series set. Furthermore, H i j is an element of matrix H
Fig. 2. Snapshots for some realizations of the dynamics at (a) primary infection period, (b) first, (c) second, (d) third region of clinical latency and (e) AIDS. In the figures, healthy cells are presented in blue, the infected-1 in yellow, the infected-2 in green and the dead cells in red. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)
JID:PLA
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.3 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
3
where i = 1, . . . , M and j = 1, . . . , T . From the Eq. (1), we can calculate the eigenvalues λi and eigenvectors vi of cross-correlation matrix Cvi = λi vi . Moreover, the eigenvalues are ordered according to size, as λ1 λ2 λ3 · · · λ M , and their distributions are directly related to the amount of temporal series set correlations. Hence, we define a nonrandom level repulsion parameter in function of eigenvalues, which is given by
si =
λ i +1 − λ i . λi +1 − λi
(3)
The parameter si is known as nearest-neighbor spacing caused by a specific correlation structure of temporal series set and the λi +1 − λi is an average over a cross-correlation matrices ensemble. In the framework of RMT, the nearest-neighbor spacing distribution can be described by Brody distribution [17]
β
β+1
Pβ (si ) = c (β)(1 + β)si exp −c (β)si where
c (β) =
β +2 β +1
,
(4) Fig. 3. Nearest-neighbor spacing distributions of the cross-correlation matrices ensemble for the first three nearest-neighbor spacings.
β+1 ,
β = [0, 1].
(5)
The Brody distribution has two important limits, when β goes to 0 or 1. In the first case (β = 0) the Eq. (4) holds Poisson distribution P0 (si ) = e −si , which means that there is no correlation structure of temporal series set. However, for second case (β = 1) the Eq. (4) holds Gaussian orthogonal ensemble (GOE) distribution
P1 (si ) =
π 2
π
si exp −
4
s2i ,
(6)
which means that there is a maximum correlation structure of temporal series set. Furthermore, Brody distribution smoothing crossing from Poisson to GOE and can give intermediate grades of correlation structure of temporal series set at an estimation of β from 0 to 1. Beyond nearest-neighbor spacing distribution, we also can use the next nearest-neighbor spacing distribution to characterize the statistic of eigenvalues fluctuations. The next nearest-neighbor spacing is defined as following
Si =
(7)
If the nearest-neighbor spacing distribution is given by GOE the next nearest-neighbor spacing distribution is
P( S i ) =
2
36 π 3
β
β+1
Pβ (si ) = C (β)(1 + β)si exp −C (β)si
,
(9)
where
λ i +2 − λ i . λi +2 − λi
18
and d), which will be associated with the distributions of RMT and with the four regions of time dynamics evolution of HIV infection presented in Figs. 1(a, b, c and d), respectively. Firstly, with the aim of connecting the nearest-neighbor spacing statistic with RMT, we fitted the data set of histograms using the Brody distribution, Eq. (4). The Brody distribution does not fit the data set of the high order nearest-neighbors spacings showed in Figs. 4(a) and 4(c) (dashed line), while it fits the ones presented in Figs. 4(b) and 4(d) satisfactorily. Furthermore, we estimate β = 0.09 for the data in Fig. 4(b), which means the distribution is a Poisson and also that the data set of nearest-neighbor spacing lost any correlation, while we estimate β = 1.02 for the data in Fig. 4(d), which means the distribution is a GOE and that the data set of nearest-neighbor spacing has maximum correlation. As the Brody distribution cannot be used to fit the data set of Figs. 4(a) and 4(c), we have gotten a ansatz which is given by
S 4i exp −
64 9π
S 2i ,
(8)
which is known as Gaussian symplectic ensemble (GSE) in the framework of RMT.
C (β) =
β +2 β +1
(1 + β) ,
β = [0, 1].
By using the ansatz Eq. (9) to fit the data set of Figs. 4(a) and 4(c), we have obtained a great agreement (solid line). This means that the data set of Figs. 4(a) and 4(c) have a degree of correlation. The estimation of β in both cases is β ≈ 1 (0.99 and 0.95 respectively). Replacing β = 1 in Eq. (9) follows that
π
P(si ) = π si exp −
4. Results and discussion In this section, we will use the cross-correlation method described above to study the time dynamics evolution of HIV infection statistic. Using the HIV model of section 2, an ensemble composed by 3260 cross-correlation matrices C with dimension 103 × 103 was built, in which each row of the matrix corresponds to the time evolution of the density of infected cells of one individual. Let’s start plotting some nearest-neighbors spacing histograms from the cross-correlation matrices ensemble, as can be seen in the Figs. 3 and 4. The histograms of Fig. 3 for the first three nearest-neighbors spacings do not show stability, which means that they are not statistically similar. However, we can identify four stable classes of histograms with different statistics, Figs. 4(a, b, c
β+1
2
s2i ,
(10)
which is slightly different from GOE, Eq. (6). Hence, we can affirm that the nearest-neighbor spacing statistic of Figs. 4(a) and 4(c) is quite different of the ones of Fig. 4(d). Furthermore, we also built the histograms of next nearestneighbor spacing as can be seen in the Fig. 5. We have compared the data set of histograms with GSE and Poisson distributions. As expected, the results showed in Figs. 5(b) and 5(d) present a great agreement with Poisson and GSE distributions, respectively, while the ones showed in Figs. 5(a) and 5(c) were not fitted by GSE (dashed line). However, the next nearest-neighbor spacing data set of Figs. 5(a) and 5(c) can be fitted by following distribution (solid line)
P( S i ) = c¯ 1 S 4i exp −¯c 2 S 2i ,
(11)
JID:PLA
4
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.4 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
Fig. 4. Nearest-neighbor spacing distributions of the cross-correlation matrices ensemble for high order nearest-neighbors spacing. In parts (a) and (c) the data set were fitted by Eq. (9) (solid line), furthermore the dashed line shows that the Brody distribution Eq. (4) does not fit the data set. In parts (b) and (d) we used Brody distribution Eq. (4) (solid line).
Fig. 5. Next nearest-neighbor spacing distributions of the cross-correlation matrices ensemble for high order next nearest-neighbors spacing. In parts (a) and (c) the data set were fitted by Eq. (11) (solid line), furthermore the dashed line shows that GSE distribution is not in agreement with data set as it was expected. The solid line is the Poisson distribution in part (b) and the GSE distribution (Eq. (8)) in part (d).
where c¯ 1,2 are fit parameters. The results showed in Fig. 5 are in agreement with what is expected for next nearest-neighbor spacing statistic in the framework of RMT. We mention that there is not a direct association between different ranges of i and times along the evolution. However, we can do an indirect association analyzing the formation and destruction of structures in Fig. 2. In Fig. 2(a) we show a snapshot of the dynamics that represents the primary infection period (Fig. 1(PI)) which is characterized by the absence of structures. After the pri-
mary infection, the clinical latency phase is characterized by the formation of structures as shown in Fig. 2(b), corresponding to the region of Fig. 1(a). Those structures completely fill the lattice after approximately 200 weeks (Fig. 2(c)), which is illustrated by the region of Fig. 1(b). Thereafter, the structures begin to be destroyed at the end of the clinical latency phase (Fig. 2(d)), and this period is represented by the region of Fig. 1(c). Only when the structures are completely destroyed the dynamics gives rise to the onset of AIDS (Fig. 2(e)), illustrated by the region of Fig. 1(d). We can conclude
JID:PLA
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.5 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
that the formation and destruction of structures is fundamental in the Zorzenon dos Santos and Coutinho model of evolution dynamics of HIV, hence it is expected that they will affect the statistic of nearest-neighbor spacing. This expectation was hold with five different kinds of histograms of nearest-neighbor spacing that are shown in the Figs. 3 and 4(a, b, c and d). Hence, we associate them with the formation and destruction of structures in the Figs. 2(a, b, c, d and e), respectively. Summarizing, we have used the tools of the RMT to investigate the statistical behavior of infected cells during the evolution of HIV infection by using an ensemble of 3260 cross-correlation matrices. From the nearest-neighbor spacing histograms obtained from the ensemble, Figs. 4(a, b, c and d), we could identify four distinct regimes of the evolution of HIV infection (as shown in Fig. 1(a, b, c and d)). We verified that at the beginning of the so-called clinical latency phase (Fig. 1(a)) the concentration of infected cells grows slowly and evolves in a correlated way (Fig. 4(a)). This regime is followed by another one (Fig. 1(b)) in which the correlation is lost (Fig. 4(b)) and that in turn leads the system to a regime in which the increase of infected cells is faster (Fig. 1(c)) and correlated (Fig. 4(c)). In the final phase, the one in which AIDS is stablished (Fig. 1(d)), the system presents maximal correlation as demonstrated by GOE distribution (Fig. 4(d)). We verified that the dynamics of the model drastically affects the statistic of the nearest-neighbor spacing, so we expect that the statistic analysis could be used as a measure of the efficiency of antiretroviral therapies of HIV infection. The results showed in this manuscript can contribute to a better understanding of the evolution of HIV infection and future investigations can include the application of tools of RMT to dynamics of HIV infection under
5
antiretroviral therapy and also to other problems in epidemiology and population dynamics. Acknowledgements The authors acknowledge S. Coutinho for useful discussions. The authors are supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and by Fundação de Amparo a Ciência e Tecnologia do Estado de Pernambuco (FACEPE). References [1] R.M.Z. dos Santos, S. Coutinho, Phys. Rev. Lett. 87 (2001) 168102. [2] P.H. Figueirêdo, S. Coutinho, R.M.Z. dos Santos, Physica A 387 (2008) 6545. [3] G. Solovey, F. Peruani, S.P. Dawson, R.M.Z. dos Santos, Physica A 343 (2004) 543. [4] A.S. Perelson, P.W. Nelson, SIAM Rev. 41 (1999) 3. [5] Z. Hel, J.R. McGhee, J. Mestecky, Trends Immunol. 27 (6) (2006) 274. [6] O.J. Cohen, G. Pantaleo, G.K. Lam, A.S. Fauci, Springer Semin. Immunopathol. 18 (3) (1997) 305. [7] E.P. Wigner, SIAM Rev. 9 (1967) 1. [8] M.L. Mehta, Random Matrices, Academic Press, 2004. [9] R. Couillet, M. Debbah, Random Matrix Methods for Wireless Communications, Cambridge University Press, 2011. [10] O. Bohigas, M.J. Giannoni, C. Schmit, Phys. Rev. Lett. 52 (1984) 1. [11] R. Potestio, F. Caccioli, P. Vivo, Phys. Rev. Lett. 103 (2009) 268101. [12] A. Agrawal, C. Sarkar, S.K. Dwivedi, N. Dhasmana, S. Jalan, Physica A 404 (2014) 359–367. [13] F. Luo, et al., Phys. Lett. A 357 (2006) 420. [14] M. Müller, et al., Phys. Rev. E 74 (2006) 041119. [15] L. Laloux, P. Cizeau, J.P. Bouchaud, M. Potters, Phys. Rev. Lett. 83 (1999) 1467. [16] V. Plerou, P. Gopikrishnan, B. Rosenow, L.A.N. Amaral, H.E. Stanley, Phys. Rev. Lett. 83 (1999) 1471. [17] T.A. Brody, Lett. Nuovo Cimento 7 (1973) 482.
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.1 (1-5)
Physics Letters A ••• (••••) •••–•••
Contents lists available at ScienceDirect
Physics Letters A www.elsevier.com/locate/pla
Statistical behavior of time dynamics evolution of HIV infection Ramón E.R. González, Iury A.X. Santos, Marcos G.P. Nunes, Viviane M. de Oliveira ∗ , Anderson L.R. Barbosa ∗ Departamento de Física, Universidade Federal Rural de Pernambuco, 52171-900, Recife, Pernambuco, Brazil
a r t i c l e
i n f o
Article history: Received 15 May 2017 Received in revised form 14 July 2017 Accepted 17 July 2017 Available online xxxx Communicated by C.R. Doering Keywords: Random matrix theory Cellular automata HIV infection
a b s t r a c t We use the tools of the random matrix theory (RMT) to investigate the statistical behavior of the evolution of human immunodeficiency virus (HIV) infection. By means of the nearest-neighbor spacing distribution we have identified four distinct regimes of the evolution of HIV infection. We verified that at the beginning of the so-called clinical latency phase the concentration of infected cells grows slowly and evolves in a correlated way. This regime is followed by another one in which the correlation is lost and that in turn leads the system to a regime in which the increase of infected cells is faster and correlated. In the final phase, the one in which acquired immunodeficiency syndrome (AIDS) is stablished, the system presents maximum correlation as demonstrated by GOE distribution. © 2017 Elsevier B.V. All rights reserved.
1. Introduction Since the first cases of human immunodeficiency virus (HIV) infection were diagnosed many investigations have been carried out in order to understand the evolution of the infection [1–6]. The dynamics of HIV infection is characterized by three distinct phases which are related to the evolution of the degree of infectability. In the first one, the primary infection, a peak in the infection level is observed around 5 weeks after the contagion. This phase lasts from 10 to 12 weeks. After this period the immune system presents false signs of recovery. The second phase, the clinical latency period, is characterized by a very low viral concentration followed by a gradual growth which leads to the impairment of essential functions of the immune system. This period may last for 2 to 10 years. After this period the probability of the patient developing acquired immunodeficiency syndrome (AIDS) is very large. Several mathematical models have been proposed to study the dynamics of HIV infection [4]. Most of the investigations are based on the use of differential equations but models in which the spatial structure is considered have also been developed. In this sense, Zorzenon dos Santos and Coutinho [1] used a cellular automata model to study the evolution of HIV infection and the onset of AIDS. Their model propose a local interaction among the different types of cells (healthy, infected and dead), while the lymph node is modeled by using a square lattice with periodic boundary condi-
*
Corresponding authors. E-mail addresses: [email protected] (V.M. de Oliveira), [email protected] (A.L.R. Barbosa). http://dx.doi.org/10.1016/j.physleta.2017.07.022 0375-9601/© 2017 Elsevier B.V. All rights reserved.
tions and Moore neighborhood. The model was able to reproduce the time scales of the three phases of the infection. The random matrix theory (RMT) was proposed by Wigner and Dyson for the study of the statistical properties of the spectrum of complex nuclei [7,8]. This approach has been successfully applied to the study of many subjects such as wireless communications [9], chaotic systems [10], protein dynamics [11,12] and structure and interactions of biological networks [13]. Furthermore, the RMT can also be applied to multivariate time series as electroencephalographic recordings with epileptic events [14] and financial time series [15,16], which motivated us using this approach to the time dynamics of HIV infection. In this work we use the cellular automata model introduced by Zorzenon dos Santos and Coutinho [1] to study the statistic behavior of time dynamics of HIV infection by means of the RMT. In Section 2, we describe the cellular automata model used to study the evolution of HIV infection and we also present the time evolution of density of infected cells as generated by the model. A brief review about equal-time cross-correlation matrix and RMT is made in Section 3. Finally, in Section 4 we show our results and discussion. 2. Cellular automata model for the dynamics of HIV infection In the cellular automata model used by Zorzenon dos Santos and Coutinho [1] the square lattice is composed by cells that can be in one of four possible states: healthy, infected-1, infected-2 and dead. In the infected-1 state, the infected cells are able to spread the infection. The state infected-2 corresponds to infected cells which have already been recognized by the immune system.
JID:PLA
2
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.2 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
In the initial configuration, the system is composed by healthy cells and also by a small fraction of infected cells, p H I V , which are randomly distributed on the lattice. The infection process is governed by four rules: (1) A healthy cell will become an infected-1 cell if there is at least one infected-1 cell in its neighborhood (Moore neighborhood), or if there are at least R infected-2 cells in its neighborhood; (2) After spend τ time steps in the same state, an infected-1 cell will become an infected-2 cell, in which τ is the time necessary for the immune response; (3) An infected-2 cell will become a dead cell in the next time step; (4) A dead cell will become a healthy cell with probability p repl . If it happens, the healthy cell will become an infected-1 cell with probability p inf ec . After the update of the lattice according to the rules described above one time step is counted which corresponds to one week. In Fig. 1 we show the density of infected cells as a function of time, obtained from the model of Zorzenon dos Santos and Coutinho [1], in which the initial conditions are p H I V = 0.05, p inf ec = 10−5 and p repl = 0.99, τ = 4 and R = 4. We have used a lattice of size equal to L = 700 and the results are obtained from averages over 500 distinct simulations. We can clearly notice a qualitative behavior very similar to the one presented by the dynamics of the HIV infection. In the period between 0 and 12 weeks we verify the first phase of the dynamics (primary infection), that is characterized by a peak of the infected cells. In the second phase, the clinical latency, we observe three regions which are delimited in Fig. 1: In region (a) the concentration of infected
cells is very small and presents a very slow increment with time; in the second regime (b) we verify an augment of the fluctuation of the concentration of infected cells in which the curve presents a transition region with an inflection point leading to another regime (c) that is characterized by a faster increasing of the concentration of infected cells leading to the development of AIDS. In the third phase (d) we observe the onset of AIDS. We can better understand the behavior of the system by looking at the snapshots of five distinct time steps of some realizations of the dynamics. In Fig. 2(a) we can visualize the system during the primary infection period. After the introduction of the infected cells, pulses of these cells with width (τ + 1) are generated and propagate in all directions. Around the fifth week of the infection process one can notice a peak in the concentration of infected cells and almost all the sites of the lattice are composed by cells of type infected-1 and infected-2 (which are represented in yellow and green, respectively). From this time the concentration of infected cells decreases and the end of the primary infection phase is observed when the number of time steps is equal to 2(τ + 1). After this initial phase new infected cells will be originated when dead cells are replaced by infected-1 cells. They generate structures like the ones showed in Fig. 2(b), which spread over the lattice in a much longer time that the one observed in primary infection. As long as these structures increase they can occupy the whole lattice (see Figs. 2(c) and 2(d)), incrementing the number of infected cells until they reach a stationary state. At the end of this long period in which the increase of these structures is observed, we verify another pattern in which the distinct cells are randomly distributed in the lattice with the predominance of the infected cells (Fig. 2(e)), characterizing the onset of AIDS. 3. Cross-correlation and random matrix theory In this section, we will do a brief review about equal-time cross-correlation matrix [14] and RMT [8]. Their concepts will be applied to study the time dynamics evolution of HIV infection statistic in the next section. The equal-time cross-correlation matrix C is obtained from an M-dimensional temporal series set H i j with length T , where i = 1, . . . , M and j = 1, . . . , T with M ≤ T . The correlation matrix is given by
C=
1 T
H Ht ,
(1)
where index t is a transposed of matrix and the H matrix elements are built from temporal series set as following
Fig. 1. Time evolution of the density of infected cells obtained from the simulations. The results are obtained from averages over 500 simulations. We used lattices with size L = 700, p H I V = 0.05, p inf ec = 10−5 , p repl = 0.99, τ = 4 and R = 4. PI means primary infection.
Hij =
Hij − Hi
σHi
.
(2)
The H i and σ H i are the average and standard deviation of temporal series set. Furthermore, H i j is an element of matrix H
Fig. 2. Snapshots for some realizations of the dynamics at (a) primary infection period, (b) first, (c) second, (d) third region of clinical latency and (e) AIDS. In the figures, healthy cells are presented in blue, the infected-1 in yellow, the infected-2 in green and the dead cells in red. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)
JID:PLA
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.3 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
3
where i = 1, . . . , M and j = 1, . . . , T . From the Eq. (1), we can calculate the eigenvalues λi and eigenvectors vi of cross-correlation matrix Cvi = λi vi . Moreover, the eigenvalues are ordered according to size, as λ1 λ2 λ3 · · · λ M , and their distributions are directly related to the amount of temporal series set correlations. Hence, we define a nonrandom level repulsion parameter in function of eigenvalues, which is given by
si =
λ i +1 − λ i . λi +1 − λi
(3)
The parameter si is known as nearest-neighbor spacing caused by a specific correlation structure of temporal series set and the λi +1 − λi is an average over a cross-correlation matrices ensemble. In the framework of RMT, the nearest-neighbor spacing distribution can be described by Brody distribution [17]
β
β+1
Pβ (si ) = c (β)(1 + β)si exp −c (β)si where
c (β) =
β +2 β +1
,
(4) Fig. 3. Nearest-neighbor spacing distributions of the cross-correlation matrices ensemble for the first three nearest-neighbor spacings.
β+1 ,
β = [0, 1].
(5)
The Brody distribution has two important limits, when β goes to 0 or 1. In the first case (β = 0) the Eq. (4) holds Poisson distribution P0 (si ) = e −si , which means that there is no correlation structure of temporal series set. However, for second case (β = 1) the Eq. (4) holds Gaussian orthogonal ensemble (GOE) distribution
P1 (si ) =
π 2
π
si exp −
4
s2i ,
(6)
which means that there is a maximum correlation structure of temporal series set. Furthermore, Brody distribution smoothing crossing from Poisson to GOE and can give intermediate grades of correlation structure of temporal series set at an estimation of β from 0 to 1. Beyond nearest-neighbor spacing distribution, we also can use the next nearest-neighbor spacing distribution to characterize the statistic of eigenvalues fluctuations. The next nearest-neighbor spacing is defined as following
Si =
(7)
If the nearest-neighbor spacing distribution is given by GOE the next nearest-neighbor spacing distribution is
P( S i ) =
2
36 π 3
β
β+1
Pβ (si ) = C (β)(1 + β)si exp −C (β)si
,
(9)
where
λ i +2 − λ i . λi +2 − λi
18
and d), which will be associated with the distributions of RMT and with the four regions of time dynamics evolution of HIV infection presented in Figs. 1(a, b, c and d), respectively. Firstly, with the aim of connecting the nearest-neighbor spacing statistic with RMT, we fitted the data set of histograms using the Brody distribution, Eq. (4). The Brody distribution does not fit the data set of the high order nearest-neighbors spacings showed in Figs. 4(a) and 4(c) (dashed line), while it fits the ones presented in Figs. 4(b) and 4(d) satisfactorily. Furthermore, we estimate β = 0.09 for the data in Fig. 4(b), which means the distribution is a Poisson and also that the data set of nearest-neighbor spacing lost any correlation, while we estimate β = 1.02 for the data in Fig. 4(d), which means the distribution is a GOE and that the data set of nearest-neighbor spacing has maximum correlation. As the Brody distribution cannot be used to fit the data set of Figs. 4(a) and 4(c), we have gotten a ansatz which is given by
S 4i exp −
64 9π
S 2i ,
(8)
which is known as Gaussian symplectic ensemble (GSE) in the framework of RMT.
C (β) =
β +2 β +1
(1 + β) ,
β = [0, 1].
By using the ansatz Eq. (9) to fit the data set of Figs. 4(a) and 4(c), we have obtained a great agreement (solid line). This means that the data set of Figs. 4(a) and 4(c) have a degree of correlation. The estimation of β in both cases is β ≈ 1 (0.99 and 0.95 respectively). Replacing β = 1 in Eq. (9) follows that
π
P(si ) = π si exp −
4. Results and discussion In this section, we will use the cross-correlation method described above to study the time dynamics evolution of HIV infection statistic. Using the HIV model of section 2, an ensemble composed by 3260 cross-correlation matrices C with dimension 103 × 103 was built, in which each row of the matrix corresponds to the time evolution of the density of infected cells of one individual. Let’s start plotting some nearest-neighbors spacing histograms from the cross-correlation matrices ensemble, as can be seen in the Figs. 3 and 4. The histograms of Fig. 3 for the first three nearest-neighbors spacings do not show stability, which means that they are not statistically similar. However, we can identify four stable classes of histograms with different statistics, Figs. 4(a, b, c
β+1
2
s2i ,
(10)
which is slightly different from GOE, Eq. (6). Hence, we can affirm that the nearest-neighbor spacing statistic of Figs. 4(a) and 4(c) is quite different of the ones of Fig. 4(d). Furthermore, we also built the histograms of next nearestneighbor spacing as can be seen in the Fig. 5. We have compared the data set of histograms with GSE and Poisson distributions. As expected, the results showed in Figs. 5(b) and 5(d) present a great agreement with Poisson and GSE distributions, respectively, while the ones showed in Figs. 5(a) and 5(c) were not fitted by GSE (dashed line). However, the next nearest-neighbor spacing data set of Figs. 5(a) and 5(c) can be fitted by following distribution (solid line)
P( S i ) = c¯ 1 S 4i exp −¯c 2 S 2i ,
(11)
JID:PLA
4
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.4 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
Fig. 4. Nearest-neighbor spacing distributions of the cross-correlation matrices ensemble for high order nearest-neighbors spacing. In parts (a) and (c) the data set were fitted by Eq. (9) (solid line), furthermore the dashed line shows that the Brody distribution Eq. (4) does not fit the data set. In parts (b) and (d) we used Brody distribution Eq. (4) (solid line).
Fig. 5. Next nearest-neighbor spacing distributions of the cross-correlation matrices ensemble for high order next nearest-neighbors spacing. In parts (a) and (c) the data set were fitted by Eq. (11) (solid line), furthermore the dashed line shows that GSE distribution is not in agreement with data set as it was expected. The solid line is the Poisson distribution in part (b) and the GSE distribution (Eq. (8)) in part (d).
where c¯ 1,2 are fit parameters. The results showed in Fig. 5 are in agreement with what is expected for next nearest-neighbor spacing statistic in the framework of RMT. We mention that there is not a direct association between different ranges of i and times along the evolution. However, we can do an indirect association analyzing the formation and destruction of structures in Fig. 2. In Fig. 2(a) we show a snapshot of the dynamics that represents the primary infection period (Fig. 1(PI)) which is characterized by the absence of structures. After the pri-
mary infection, the clinical latency phase is characterized by the formation of structures as shown in Fig. 2(b), corresponding to the region of Fig. 1(a). Those structures completely fill the lattice after approximately 200 weeks (Fig. 2(c)), which is illustrated by the region of Fig. 1(b). Thereafter, the structures begin to be destroyed at the end of the clinical latency phase (Fig. 2(d)), and this period is represented by the region of Fig. 1(c). Only when the structures are completely destroyed the dynamics gives rise to the onset of AIDS (Fig. 2(e)), illustrated by the region of Fig. 1(d). We can conclude
JID:PLA
AID:24627 /SCO Doctopic: Statistical physics
[m5G; v1.221; Prn:27/07/2017; 9:52] P.5 (1-5)
R.E.R. González et al. / Physics Letters A ••• (••••) •••–•••
that the formation and destruction of structures is fundamental in the Zorzenon dos Santos and Coutinho model of evolution dynamics of HIV, hence it is expected that they will affect the statistic of nearest-neighbor spacing. This expectation was hold with five different kinds of histograms of nearest-neighbor spacing that are shown in the Figs. 3 and 4(a, b, c and d). Hence, we associate them with the formation and destruction of structures in the Figs. 2(a, b, c, d and e), respectively. Summarizing, we have used the tools of the RMT to investigate the statistical behavior of infected cells during the evolution of HIV infection by using an ensemble of 3260 cross-correlation matrices. From the nearest-neighbor spacing histograms obtained from the ensemble, Figs. 4(a, b, c and d), we could identify four distinct regimes of the evolution of HIV infection (as shown in Fig. 1(a, b, c and d)). We verified that at the beginning of the so-called clinical latency phase (Fig. 1(a)) the concentration of infected cells grows slowly and evolves in a correlated way (Fig. 4(a)). This regime is followed by another one (Fig. 1(b)) in which the correlation is lost (Fig. 4(b)) and that in turn leads the system to a regime in which the increase of infected cells is faster (Fig. 1(c)) and correlated (Fig. 4(c)). In the final phase, the one in which AIDS is stablished (Fig. 1(d)), the system presents maximal correlation as demonstrated by GOE distribution (Fig. 4(d)). We verified that the dynamics of the model drastically affects the statistic of the nearest-neighbor spacing, so we expect that the statistic analysis could be used as a measure of the efficiency of antiretroviral therapies of HIV infection. The results showed in this manuscript can contribute to a better understanding of the evolution of HIV infection and future investigations can include the application of tools of RMT to dynamics of HIV infection under
5
antiretroviral therapy and also to other problems in epidemiology and population dynamics. Acknowledgements The authors acknowledge S. Coutinho for useful discussions. The authors are supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and by Fundação de Amparo a Ciência e Tecnologia do Estado de Pernambuco (FACEPE). References [1] R.M.Z. dos Santos, S. Coutinho, Phys. Rev. Lett. 87 (2001) 168102. [2] P.H. Figueirêdo, S. Coutinho, R.M.Z. dos Santos, Physica A 387 (2008) 6545. [3] G. Solovey, F. Peruani, S.P. Dawson, R.M.Z. dos Santos, Physica A 343 (2004) 543. [4] A.S. Perelson, P.W. Nelson, SIAM Rev. 41 (1999) 3. [5] Z. Hel, J.R. McGhee, J. Mestecky, Trends Immunol. 27 (6) (2006) 274. [6] O.J. Cohen, G. Pantaleo, G.K. Lam, A.S. Fauci, Springer Semin. Immunopathol. 18 (3) (1997) 305. [7] E.P. Wigner, SIAM Rev. 9 (1967) 1. [8] M.L. Mehta, Random Matrices, Academic Press, 2004. [9] R. Couillet, M. Debbah, Random Matrix Methods for Wireless Communications, Cambridge University Press, 2011. [10] O. Bohigas, M.J. Giannoni, C. Schmit, Phys. Rev. Lett. 52 (1984) 1. [11] R. Potestio, F. Caccioli, P. Vivo, Phys. Rev. Lett. 103 (2009) 268101. [12] A. Agrawal, C. Sarkar, S.K. Dwivedi, N. Dhasmana, S. Jalan, Physica A 404 (2014) 359–367. [13] F. Luo, et al., Phys. Lett. A 357 (2006) 420. [14] M. Müller, et al., Phys. Rev. E 74 (2006) 041119. [15] L. Laloux, P. Cizeau, J.P. Bouchaud, M. Potters, Phys. Rev. Lett. 83 (1999) 1467. [16] V. Plerou, P. Gopikrishnan, B. Rosenow, L.A.N. Amaral, H.E. Stanley, Phys. Rev. Lett. 83 (1999) 1471. [17] T.A. Brody, Lett. Nuovo Cimento 7 (1973) 482.