Turing Patterns of Non-linear S-I Model on Random and Real-Structure Networks with Diarrhea Data

Most developed models for solving problems in epidemiology use deterministic approach. To cover the lack of spatial sense in the method, one uses statistical modeling, reaction-diffusion in continuous medium, or multi-patch model to depict epidemic activities in several connected locations. Here, we show that an epidemic model that is set as an organized system on networks can yield Turing patterns and other interesting behaviors that are sensitive to the initial conditions. The formed patterns can be used to determine the epidemic arrival time, its first peak occurrence and the peak duration. These epidemic quantities are beneficial to identify contribution of a disease source node to the others. Using a real structure network, the system also exhibits a comparable disease spread pattern of Diarrhea in Jakarta.

idea was inspired by Fisher's equation, which explained the dispersion of organism using the logistic equation 21 . Cruickshank 13 then adopted the logistic equation as the recruitment rate for the susceptible population in the S-I model. However, the dispersion of the organism only occurs for the contagious individual. It was then extended by considering the diffusion on the susceptible individuals as well 14,15 , cross-diffusion events (that involve the motion of the susceptible individuals because of the presence of contagious ones) 16,17,22 and the basic S-I-R model of Kermack and McKendrick 1 in a reaction-diffusion system 23 . The latter type of extensions shows the different kinds of spatial patterns in the spread of diseases 24 with transitions as an emergent property that can serve as a potential trend indicator 25 .
Following the inception of network studies that were inspired by the large structure of internet connections [26][27][28][29] , epidemic models with network connection to incorporate the spatial dependence were then introduced [28][29][30][31] . In this case, each node on the network represents an individual in the population and its state can change from healthy into infected depending on the connection to infected individuals with some transmission rate 30 . An epidemic model on networks using three compartments was then developed, where in this case each node represents a group of population containing susceptible, infected, and recovered individuals 31 . The model could explain the phenomenon that epidemic cases depend on social topology and may rise when an epidemic attacks individuals with high connectivities 31 . How human behavioural and social dynamics affects disease spreading and prevention in well-mixed and networked populations have also been widely studied 32 .
A different perspective was introduced by considering each node as an area and disease transmission through the population occurs in that area 28,29 . Connections between nodes, i.e. edges, represent the motion of people between the areas. This model was developed to capture disease transmission phenomena in metapopulation events 28 , that can explain the influence of connections through their coupling strength 29 . It was shown that disease transmissions on Winnipeg city and its smaller neighborhood cities could be driven by a small city with low coupling strength. Modelling in metapopulation or networks can give an understanding the effects of people's motions in disease spreading. However, recognizing disease spread pattern as the mobility consequence is challenging since it can involve many factors such as motion trajectories and directions between some locations, pattern of people movement which can represent commuter or migration activity and how fast people move from one area to another. Recently, the study on pattern formation through reaction-diffusion model arose in network organized systems 33 . The network organized reaction-diffusion model exhibits patterns which represent movement structurally easier to understand than continuous model. We may be difficult to identify the behavior of continuously area since the original reaction-diffusion model represents population density on every point. In fact, people are occupying places discretely and continuous model may be not appropriate in disease spread case. Developing model of a big region into some points offers model relevance to real problem and clearer interpretation. Therefore, exploring disease spread using reaction-diffusion model in network-organized structure is considerably more appropriate.
In this work, we investigate the possibilities of emerging patterns from the reaction-diffusion of an S-I model in networks. We consider a two compartments epidemic model with logistic growth and a not constant total population. Numerical simulations have been done to reveal the complex patterns in the S-I model that depends on diffusion coefficient parameters.
Analysis of pattern emergence presented in this paper uses a mean-field approximation, visualizing the stationary solution in the average sense [33][34][35] . Even though the diffusion coefficient is set below threshold, pattern formation can occur and will be shown through hysteresis analysis. It emerges because of an initial condition selection. By varying the initial condition of certain nodes and set the other nodes near a homogeneous equilibrium, we present simulations of differentiated nodes along with the analysis to compute their values at steady states.
Our model represents a situation where locations with similar characteristics show a variation in their disease spread patterns. The important findings of the developed model are epidemic arrival time pattern of an emerging infected node and the occurrence of first epidemic peak and its duration. Node degrees and mobility rate drive the disease to distribute broadly. We present epidemic diarrhea data of Jakarta and show that the pattern is qualitatively comparable with reaction-diffusion model.
The outline of the paper is given as follows. We start with introducing the model. Next, we discuss the system's equilibria and study their stability, which provide preceding information to pattern formation. Numerical simulations for moderate and large random networks are then performed. The S-I model shows the occurrence of hysteresis when the ratio of the diffusion coefficients between the populations is set below a critical value. A single-differentiated-node analysis is performed to understand the numerical observation. Finally, a real network structure representing districts of Jakarta is considered and epidemic diarrhea data shows a comparable pattern with the numerical results.

epidemic Model in Networks
Full model. The reaction-diffusion systems in networks 33,36,37 are described generally by where u i (t) and v i (t) are local densities of the activator and inhibitor species at node i and time t. The kinetic functions on the i th node are given by f(u i , v i ) and g(u i , v i ). The diffusion process in the system is going through the network, whose topology is defined by N number of nodes and represented by N × N adjacency matrix A with the www.nature.com/scientificreports www.nature.com/scientificreports/ element A ij equal to 1 if there is a link between i th node and j th node and 0 if there is no link between them. The degree of node i is given by . The Laplacian matrix L whose element is given by L ij = A ij − k i δ ij where δ ij is Kronecker's delta that is equal to 1 when i = j and 0 when i ≠ j. D u and D v act as the diffusion coefficients of species u and v, respectively. The coefficients D u and D v are then scaled to D v /D u = σ and D u = ε.
As a particularly interesting case of (1), here we consider a network of the nonlinear S-I model S i (t) and I i (t) are local densities of susceptible and infected individuals, respectively, in i th node at time t. The susceptible and the infected individuals are analogous to the activator and inhibitor species in system (1), where population increment happens in the susceptible compartment as disease can inhibit the growth of the population. Parameter r is the growth rate of the susceptible population. Parameters K, β, and γ are the carrying capacity, the transmission rate, and the mortality rate caused by the disease, respectively. The form of the kinetic functions f and g is due to Cruickshank et al. 13 and Sun et al. 22 with the logistic-like equation as the recruitment rate 13,14,17,22 , which usually is assumed to be constant. In this case, the transmission process occurs between susceptible and a proportion of infected, i.e. not the total, population. The interaction that occurs in the model represents a disease transmission between locations. Herein, the networks are assumed to be scale-free.
Instability of the uniform equilibria. Without diffusion in (2), i.e. ε = 0, one can easily identify that the remaining system has two equilibria. They are the disease free equilibrium E 1 = (S * , I * ) with S * = K and I * = 0, and the endemic equilibrium provided that to be biologically relevant (as outside the interval, I * < 0). Starting from the neighborhood of E 1 , the system will allow the susceptible to increase its population until it reaches the carrying capacity, i.e. the equilibrium is unstable. On the other hand, the endemic equilibrium E 2 preserves the disease to stay and infects the population for a long time, i.e. the equilibrium is stable (see Methods section). Now we consider the S-I model (2) with diffusion. We consider a scale-free network with a size of N = 300 and mean degree 〈k〉 = 10. The parameters are set to r = 0.27, K = 1000, β = 0.5, and γ = 0.25. The distribution of the node degrees of the network is shown in Fig. 1. We sort the nodes in the order of their degrees where index 1 refers to the node with the largest node degree and index N to that with the smallest node degree.
To exhibit patterns via Turing bifurcation 38 , the critical eigenvalues of the linearisation operator around an equilibrium of the model must be complex valued with positive real part. The diffusion process can only change the stability of the endemic equilibrium E 2 , not the disease free E 1 (see Methods section). Hence, below we only www.nature.com/scientificreports www.nature.com/scientificreports/ consider the former equilibrium. From the linearized equations, we obtain that to obtain unstable eigenvalues, the parameters need to satisfy the condition (β 2 − rβ − γ 2 )/β > 0.
The eigenvalues are depicted in Fig. 2 corresponding to three different values of σ. From the plots, we can see that Turing bifurcation shall occur when σ is above the critical value σ c = 2.77778. In this case, patterns will emerge in the system. However, for a fixed number of nodes N, the coupling strength ε, which is the diffusion coefficient of susceptible compartment, will play an important role. From Fig. 2b,c, Turing instability does not occur even though parameter σ has been set larger than the critical value. The coupling strength between the nodes ε cannot be too small nor too large to obtain the partition showing Turing instability. Larger size of network is necessary to yield Turing instability in Fig. 2b. While that in Fig. 2c shows declining real part of the eigenvalues. We can conclude that we may not always obtain Turing instability although the critical parameter has been met.
The critical value σ c can be obtained analytically from the characteristics polynomial equation The critical parameters σ c and θ c correspond to the condition λ(θ c ) = λ′(θ c ) = 0 and λ′′(θ c ) < 0 where prime is the derivation respects to θ. Thus, the critical parameters are given by The expression will yield σ c = 2.77778 for r = 0.27, K = 1000, β = 0.5, and γ = 0.25, which is in agreement with the numerical observation.
pattern emergence in the s-I model. For the network organization shown in Fig. 1, even though we set the parameters to be the same for all the nodes, the diffusion process can make differences for the stationary solution of system (2). Using the endemic equilibrium perturbed with a small deviation as the initial condition for the simulation, we show in Fig. 3a,b the stationary solution of system (2) using σ = 3 and in Fig. 3c,d using σ = 5. Magenta and red nodes depict the stationary solution of susceptible and infected compartments, respectively. As σ increases, the stationary solution will generally split into two clusters. However, the nodes with high connectivity will maintain its position near the weighted average of the stationary solution of each compartments. www.nature.com/scientificreports www.nature.com/scientificreports/ The behavior of the stationary solution can be explained by the mean-field approximation [33][34][35] , which is the averaged value of all the nodes. Defining the weighted average of the stationary solution as , as the global mean-field, the mean-field equation is then described by The solutions of the mean-field Eq. (7) are depicted in Fig. 3. Good agreement is clearly seen. The mean-field equations admit two solutions, which are stable and unstable. The stable and the unstable solutions are represented by dark blue line and by the light blue line, respectively.
Hysteresis of patterns. Stationary solutions of the main Eq. (2) depend on the initial condition. We will show it from a hysteresis phenomenon.
Define the amplitude A as the distance of fields in each node from the equilibrium of (2) without diffusion, i.e. Fig. 4a. The presence of a jump from the endemic equilibrium A = 0 to non-zero A as we increased σ is due to Turing instability we discussed earlier, where the endemic equilibrium changes its stability. This occurs at the critical value σ c . Normally, Turing patterns emerge as a pitchfork bifurcation. In the supercritical case, one would expect the nodes to relax to the endemic equilibrium again when parameter σ is set below the critical value. However, from Fig. 4a, we observe a hysteresis where this is not the case. Figure 4b shows a non-uniform stationary solution of (2) at σ = 2.7, which is below the critical value.  www.nature.com/scientificreports www.nature.com/scientificreports/ From the figure, we obtain that nodes with high connectivity maintain its position near the equilibrium.
Notably we see the presence of isolated nodes with a value that is relatively distinguishable and far from the average value. Such a node is referred to as a single differentiated node (SDN) 36 .
To analyze SDN, we shall assume that the nodes are mostly in the endemic equilibrium except the i th one. Then, system (2) can be reduced into where k i is the node degree. Stationary solutions of (8) are depicted in Fig. 5 for two different node degrees k i = 12 in Fig. 5a,b and k i = 5 in Fig. 5c,d. The vertical axis is σ, while the horizontal one is the dependent variables. In the figures, blue line and red dashed line represent stable and unstable solutions, respectively. The plots show that to obtain an SDN, the higher the node degree, the higher the value of σ is required. This can explain the persistence of nodes with high connectivity to maintain their position near the endemic equilibrium, while nodes with lower connectivity can vary quite far from the endemic equilibrium.
To show that, we have performed several simulations using networks with moderate size (N = 30) and carefully selected initial conditions. By initially choosing the predicted stable SDN for the i th node and endemic equilibrium for the others, we show the results in Fig. 6, where one can clearly observe the presence of SDN in the network. This hence confirms the observation of non-vanishing amplitude in Fig. 4 when the diffusion ratio σ is even below the critical value σ c . (2) is represented by the parameters ε and σ, which show the movement rate of healthy and infected populations, respectively. The disease spread is affected by the arrival of symptomatic or asymptomatic ill persons at a new destination 39 . While in previous sections, we are interested in the characteristics of stationary patterns, it is important to also study transient dynamics of the system leading to stationary states. For the sake of clarity, we now consider moderate size network. Nodes are ordered based on the size of their connectivity degree. The first and the last index indicate node with the largest and smallest node degree, respectively. Taking different values of σ, we show the resulting stationary patterns in Fig. 7. We obtain that the distribution of infected people tends to spread evenly for a smaller σ, i.e., low movement rate of infected compartment (Fig. 7a). A larger σ (Fig. 7b) makes the infected population concentrate at nodes with moderate and large degree connectivities. We have investigated system (2) using a moderate size network when initially an infection is introduced at one node. Figure 8 shows the dynamics of the system as the solution evolves in time. It is interesting to note that before reaching a stationary state, the system oscillates.

Effects of parameters on non-uniform states. Mobility in system
It is important to note that the resulting pattern is rather independent of the growth, transmission, and death rates. The carrying capacity does not influence the longterm pattern. It only affects the magnitude of the compartments. Nevertheless, those parameters has a significant impact on the so-called first epidemic peak which will be explained in the next sub-section. Considering a larger K yields a slower first epidemic peak as illustrated in Fig. 8b,d. This finding shows the linearity between carrying capacity and first peak occurrence, i.e., the disease needs more time to infect large population. epidemic arrival and peak times. The epidemic arrival time can be defined as the time when an infected person is detected on a disease free node. This event can be identified when there is only one node which has been infected while the other nodes are free from the disease and the infected node spreads the disease to the other nodes after that. Source of infection node becomes important in distributing the disease among the nodes. On the average, the disease is quickly distributed when node with large connectivity is infected initially. Degree of node in this study has an equivalent role to the effective distance in Brockmann and Helbing article 40 . In that article, an effective distance can be used www.nature.com/scientificreports www.nature.com/scientificreports/ to estimate the epidemic arrival time in new locations. Here, the movement rates between nodes are assumed to be equal and the significance of mobility in the model is explained through node connectivity. Besides connection factor, the mobility of an infected person also has a role in accelerating the transmission to a new location (Fig. 8e). A larger movement rate causes nodes with larger connectivity to have a similar arrival time (Fig. 8e). As discussed by Arino and Portet 29 , nodes with high connectivity (small distance with high movements) can be indistinguishable; in our case, the distance is represented by the node degree. It is noted that restraining a sick person to travel to another location becomes important when the disease has a potential to be transmitted easily and globally 39 .
The first epidemic peak is defined as the period when the size of infected population is about 10% of total population in the early time of first infection and back to the same size after a certain duration. Population size shows a little significant impact on the arrival time but it has a connection with peak occurrence. First epidemic peak occurs as the disease emerges on that location. Larger population size requires a longer time to reach its peak. On the other hand, a larger mobility rate delays the peak and reduces its duration to occur although it can increase the size of infected population. Table 1 explains the overall average and deviation for epidemic time of arrival, first epidemic peak occurrence and duration.
Real-structure network patterns and exploration using real data. In this paper, we also consider a real network structure based on a map. Here, we have chosen Jakarta as our study area to explore disease occurrence. We create a network from the map by connecting every two districts that are adjacent to each other. Practically we examine that condition by finding an edge to overlap between the two districts. The connection between districts in Jakarta is depicted in Fig. 9a where we obtain a network with 42 nodes.
Our real-structure network has different distribution of node connectivity from scale-free networks considered. Figure 9b depicts the degree of node in Jakarta network that are less vary than scale-free network. Each node node in our real-structure network has at least two and no more than eight connections. On the other hand, we can find some nodes with relatively very large connectivity in scale-free network. Our network generally has similar simulation results compared to scale-free network. Although the epidemic time of arrival has an increasing trend, nodes with larger connectivity do not necessarily accelerate the emergence of disease in new locations (Fig. 10a).
Using data obtained from surveilans-dinkesdki.net 41 , we compare a real epidemic condition in Jakarta with theory in this paper. Diarrhea is an example which disease transmission between several locations can be www.nature.com/scientificreports www.nature.com/scientificreports/ considered as reaction-diffusion process on network-organized system. Figure 11a shows diarrhea cases in 2015 for every districts. The pattern shows that the disease occurs in a year. We simulate model (2) using parameters in Fig. 2 with carrying capacity K = 242332. The population in Jakarta in 2015 was around 10 million based on  www.nature.com/scientificreports www.nature.com/scientificreports/ jakarta.bps.go.id 42 . Thus, the total population is divided by 42 districts and yielding around 242,332. We assume the carrying capacity for each districts is the same to reduce the complexity. We use σ = 3 and σ = 5 to compare the resulted patterns to diarrhea pattern in Jakarta. Stationary patterns from the model with σ = 5 in Fig. 11b,c exhibits a comparable pattern to the diarrhea pattern in Jakarta even though the magnitude of simulation result cannot be compared. Initial conditions for the model are selected from perturbed endemic equilibrium of system (3) and one infected in one node while the other nodes are disease free. www.nature.com/scientificreports www.nature.com/scientificreports/ Nodes which indicate high cases of diarrhea (average value larger than third quartile of 42 nodes average cases) are 3rd, 9th, 10th, 13th, 24th, 25th, 29th, 33rd, 34th and 42nd node. Simulations with 43 various of initial conditions have been explored to find nodes which indicate high infection. All nodes are classified into four class to see the potential to have high cases of disease. Class of high, medium-high, medium-low and low indicate nodes which occur very often (more than 16 times), often (between 8 and 16 times), rarely (between 3 and 8 times) and very rarely (below 3 times) with high cases in simulation, respectively. The classification result is showed in Table 2. Although the network model shows different nodes of high infection compared with real data, we can learn from this model that nodes with connection number near to average degree of our real-structure network are more likely to show high infection. Nodes with lower or higher degree exhibit to be more difficult to classified as nodes with high infection. More adjustment on the network, such as includes weight on the adjacency matrix to represent real situation on each node, is necessary to perform better result to compare with real data.

Conclusion
Studying disease transmission using reaction-diffusion in network-organized system shows different perspective in understanding epidemic patterns through deterministic mathematical model. Well-known differential equations to study disease spread is limited to temporal sense explanation on one site. By the time goes on, the epidemic dynamics cannot be restricted to a certain place. The incidence which occurs in one site propagates to different sites. One possibility which can create this condition is human mobility. Further studies in disease spreads reveal the involvement of human mobility in broad infection 28,29,39,43,44 , the mobility can be considered as a diffusion-like process in a fixed structure network.
Our results show that the endemic equilibrium can yield patterns by tuning the mobility parameter of the infected compartment above a critical value. However, Turing instability cannot emerge when mobility for susceptible compartment is wrongly selected. The stationary behavior of the system on the long run can be estimated by studying a mean-field approximation. On the other hand, time independent solution can be different from the mean-field sense. Changes of pattern from the average are shown by hysteresis and SDN events by selecting certain initial conditions. Mobility rate σ and carrying capacity K determine the infected people distribution among the nodes. Bigger value of the rate creates the infected population to be accumulated on several nodes.    www.nature.com/scientificreports www.nature.com/scientificreports/ Meanwhile, carrying capacity only affects number of victims, but different magnitudes show similar patterns. Extending the simulations by applying only one node as the infection source shows disease propagation time and its period. The arrival time is identified once the disease reaches in a node. Disease early occurrence is strongly affected by mobility rate σ and degree of nodes. Higher people movements will quicken the spread while coupling number can show the impact of source node to the other nodes. After the first occurrence propagates, it will create its first peak which be more affected by carrying capacity.
The structure of network has important role in revealing the epidemic and its behavior 40,43,44 . Networks established from real-structure arrangement of Jakarta map also is considered in this study. In general, the network does not show particular different in pattern formation, but cannot exhibit similar characteristic on epidemic time to scale-free network. Different properties of network might create distinct phenomena 40 . However, applying real-structure network on system (2) shows comparable pattern to real data by tuning parameter σ properly. This finding creates bigger opportunities to compare reaction-diffusion model in networks to real data and learn disease pattern more widely.
As in this paper we showed the appearance of Turing patterns in the model numerically, it is important to analyse them analytically. It can be done by deriving an amplitude equation using a multiple scale expansion method [45][46][47] , that can be used to reveal various patterns in the model systematically 45 . This is addressed for future work. Additionally, our work on the effects of directed network topology, that has been shown to significantly change the rule of instability occurrence 37 , in the S-I model will be reported elsewhere. that yields the polynomial characteristics λ 2 + (γ + r − β)λ + (β − γ)(γ + r − β)/β = 0. Through the Routh-Hurwitz stability criterion, the endemic equilibrium is stable when (γ + r − β) > 0 and (β − γ) (γ + r − β)/β > 0. This implies that for stability parameter β needs to be in the interval γ < β < γ + r. where λ α denotes eigenvalues of the linearized system (2). eigenvalues of the disease free equilibrium. When the Jacobian matrix of system (3) is evaluated at the disease free equilibrium, we will obtain f S = −r, f I = −β, g S = 0, and (β − γ). To obtain stable endemic equilibrium, parameter β − γ should be negative. The eigenvalues of system (11) for each Λ α are given by λ εΛ = − + α α r (1) and λ β γ σεΛ = − + α α (2) . The eigenvalues of the Laplacian matrix will be non-positive, i.e. Λ α ≤ 0. One easily can see that eigenvalues of system (11) evaluated at the disease free equilibrium will be negative. Therefore, system (2) cannot admit pattern formations around the disease free equilibrium.