The peculiar nature of Covid-19 disease is the special form of pneumonia it causes and its fast spreading in the entire population. This dreaded Covid-19 disease is caused by the SARS-CoV-2 virus, and the disease outbreak was announced as a pandemic of special health attention by the World Health Organization (WHO) on 11 March, 2020. As of June 2, 2020 (10:06am CEST), there is a total Covid-19 confirmed case of 6,140,934 with 373,548 deaths all over the world [1] . India is one of the countries which is now seriously affected by this pandemic. As of June 2, 2020 (7:45pm IST), there are a total of 2,01,341 confirmed cases with 99,135 active cases and 5,632 deaths in India [2] . The present situation of the outbreak is alarming since there is no vaccine/drug so far to cure this disease. The Indian Government, both at the central and state levels, has taken up special measures such as quarantine, social distancing and lockdown to prevent/intervene this pandemic across the country.
Mathematical epidemic models play important roles to study and predict disease dynamics as well as to implement necessary intervention strategies to control disease outbreak. Classic compartmental models such as the SI, SIS, SIR [3] , SEIR and their derived/extended models [4] have long been successfully used to study various disease transmission dynamics for different viruses such as H1N1 virus [5] , Ebola virus [6] , SARS-CoV [7] , MERS-CoV [8] , etc. With special reference to the ongoing Covid-19 pandemic, there have been attempts using statistical methods, deterministic compartmental modeling, large scale simulation to study the Covid-19 disease dynamics in order to assist in mitigating the disease outbreak in several countries across the world. There have also been various studies in India regarding this Covid-19 outbreak in the country using various mathematical models which are mostly deterministic models [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] . On the other hand, to capture the qualitative as well as quantitative real dynamical situations to intervene disease outbreak, stochastic approach needs to be employed. With the increasing capacity of modern computers, stochastic methods are gaining popularity because of being a powerful means to study and predict any complex dynamical system inherent with various environmental fluctuations/noise. In this work, we study the Covid-19 disease dynamics in India and some of its states using the classic SEQIR model from stochastic dynamics approach. Using stochastic numerical simulation, we illustrate in a very simplified manner the important impacts of social distancing and quarantine on the disease spreading mechanism specially in India and few states. We also highlight the importance of demographic stochasticity in the disease spreading dynamics.
The classical SEIR (Susceptible, Exposed, Infected and Recovered) model [4] was extended with another sensitive compartment called quarantine and the model is known as SEQIR model [16] as shown in Fig. 1 . In this model, the total population is sub-divided into five sub-populations (compartments) such as susceptible (S); exposed and infected but not detected by testing (E); self and/or institutional quarantined (Q); confirmed/reported/hospitalized infected population (I) and disease-recovered as well as population living in secured zone not affected by the Covid-19 outbreak (R). We assume uniform mixing or homogeneity in the large population and also assume the total population at any instant of time is, N = S(t) + E(t) + Q(t) + I(t) + R(t). Further, the change in population is in discrete integer amounts and the process of changing is a stochastic Markov process. Hence, the time evolution of each variable in the model should be considered in discrete and stochastic fashion [21] . In stochastic formalism of the model system [ Fig. 1 and set of coupled differential equations (1)- (5) ], the population state vector at any instant of time can be represented by, X = [S, E, Q, I, R] −1 , which undergo M = 15 reaction channels defined by, E, Q, I, R] ; i = 1, 2..., 5. Here, {a} and {b} are the sets of reactant and product molecules respectively, and {k} is the set of the classical rate constants. Further, the classical rate constants, k j can be related to the stochastic rate constants, c j by c j = k j V 1−ν , where ν is the stoichiometric ratio and V is the system size [21, 23] . This incorporates the idea of correlating fluctuations in the dynamics of the system [21] [22] [23] . Now, the reaction channels can be translated as,
where, α is the rate of disease transmission from the Susceptible to the Exposed and infected but not detected by testing class; β 1 is the rate of transition from the Susceptible to the Quarantined class; σ 1 is the rate of transition from the Susceptible to the secured zone class; r 1 is the rate of transition from the Exposed and infected but not detected by testing to the hospitalized Infected class; β 2 is the rate of transition from the Exposed and infected but not detected by testing to the Quarantine class; σ 2 is the rate of transition from Exposed and infected but not detected by testing to the secured zone class; r 2 is the rate of transition from Quarantine to hospitalized Infected class; σ 3 denotes the rate of transition from Quarantine to the secured zone class; d 2 indicates the rate of Covid-19 induced death; Λ is the rate at which new individuals enter the Indian population due to a new child-birth or immigration in the country; and d 1 denotes the natural death rate.
If the system is subjected to a certain temperature T , the trajectory of any variable in X suffers a set of random molecular events given by the set of reactions (1) and follows Brownian motion [21] [22] [23] . Further, every time any one of the reaction sets is encountered, creation and annihilation of the molecular species will take place and hence, the state vector X will get changed as a function of time. Consider the state change from state X to another state X ′ during the time interval [t, t + ∆t], then the time evolution of the configurational probability of state change P (X; t) is given by the following Master equation constructed from the detailed balance equation [24] [25] [26] , ∂P (S, E, Q, I, R; t) ∂t = α(S + 1)(E − 1) P (S + 1, E − 1, Q, I, R; t) + β 1 (S + 1) P (S + 1, E, Q − 1, I, R; t) +σ 1 (S + 1) P (S + 1, E, Q, I, R − 1; t) + r 1 (E + 1) P (S, E + 1, Q, I − 1, R; t) +β 2 (E + 1) P (S, E + 1, Q − 1, I, R; t) + σ 2 (Q + 1) P (S, E, Q + 1, I, R − 1; t) +r 2 (Q + 1) P (S, E, Q + 1, I − 1, R; t) + σ 3 (I + 1) P (S, E, Q, I + 1, R − 1; t) +d 2 (I + 1) P (S, E, Q, I + 1, R; t) + Λ P (S − 1, E, Q, I, R; t) +d 1 (S + 1) P (S + 1, E, Q, I, R; t) + d 1 (E + 1) P (S, E + 1, Q, I, R; t) +d 1 (Q + 1) P (S, E, Q + 1, I, R; t) + d 1 (I + 1) P (S, E, Q, I + 1, R; t)
The Master equation (2) for complex multivariate systems is generally difficult to solve except for simple ones. However, the Master equation of any complex system can be solved numerically using the stochastic simulation algorithm (SSA) which is discussed briefly here. SSA is generally known as Doob-Gillespie algorithm, formulated by Gillespie [21, 22] based on the theoretical foundations developed by Doob JL [27, 28] and originally proposed by Kendall [29] . It is a Monte-Carlo type of algorithm, which is a non-spatial individual based analog of the Master equation incorporating all possible interactions in the system [21] . The SSA is built on two independent processes which are random viz. firing reaction and reaction time. These independent processes are realized by the generation of two uniform random numbers r 1 and r 2 which are statistically independent. The reaction time is computed using
where a i is the i th propensity function given by a i = h i c i , where, h i is the number of possible molecular combinations of i th reaction. The j th reaction will fire when it satisfies,
We present below the simulation results of the model proposed, discuss the analysis and provide possible prediction of the Covid-19 pandemic in India and few states which are affected seriously. In population dynamics, one can express V as , V = N N/V = N D ; where, N is the total population in the geographical area, and D is the population density. Taking N as constant, we can correlate the change in V as change in the D by, V ∝ 1 D .
We numerically simulate the time-evolution of the infected population I(t) in India using SSA under different conditions. The rate constants values and initial values for the simulation are taken from [16] [30] in the increase in infected population). The curves start flattening (Gompertz-Winsor nature [31] 
, which is the signature of endemic [ Fig. 2 (b) ], and the fluctuations in the dynamics due to D playing important role in intervening the disease spreading. Hence, from these results it can be predicted that the controlling population density could lead to the endemic of Covid-19 pandemic in India. In other words, the strategy to decrease in D could be isolation of susceptible and exposed from infectious people either by quarantine them (at home or isolated place) or lockdown of socially interacting places (academic institutions, offices, festivals etc).
In Fig. 2 (c) , the sensitivity of the model is studied with respect to the parameter α, which is the disease transmission rate from S to E population Fig. 1 , at a particular value V = 0.7. As the value of α increases, the population I(t) increases sharply. This corresponds to the fact that for a certain demographic region, as the transmission rate increases, due to homogeneous mixing of population, more S population gets exposed and then infected with the disease. Further, we studied the impact of the quarantine rate, β 2 on I(t), where, the quarantine rate is given by, β 2 ∼ 1 delay in quarantine . In Fig. 2 (d) , we see that as the quarantine rate increases or delay in quarantine decreases, I(t) also starts decreasing sharply. When there is no quarantine or β 2 = 0.0, the I(t) ∼ O(10 4 ). However, when β 2 = 0.2 (quarantine in 5 days), I(t) starts flattening around ∼ O(10 3 ). This result illustrates the important effect of quarantine on I(t) indicating quarantine of population E needs to be done as quickly as possible and for longer time also for intervening the disease spreading. We know that flattening the curve can prevent the burden on hospitals and health care facilities which in turn will keep the pandemic under control.
We again perform numerical simulation of the evolution of the infected population I(t) with respect to time t in five Indian states namely, Uttar Pradesh, Delhi, Kerala, Maharashtra and West Bengal. These states are comparatively highly affected by Covid-19 pandemic and densely populated states among other Indian states. The rate constant values and initial values for the simulation of all five states are taken from [16] after verification with the present case. Fig. 2 (g) for V = 100.0, it can be seen that if the exposed E population is quarantined in ten days, the disease outbreak is relatively controlled. . We also observe that the flattening of the curve is achieved at earlier times of the disease outbreak if the values of V is increased. This could be due to the fact that the number density n of the population decreases when V is increased as the total population N is fixed. The decrease in n may be considered due to policy like social distancing, isolation etc. Our simulation study also clearly shows that social distancing plays a crucial role in early times of the epidemic for proper intervention. Again, in Fig. 3 (b) and (c), the time evolution of I(t) is studied at different quarantine rates at two different values V = 1 and V = 5 [D = N 5 ]. For the same parameter values, at V = 1 [D = N ], I(t) shows an exponential increase whereas when it is at V = 5 [D = N 5 ], I(t) already starts flattening and decreasing (also see 3(a)). In Fig. 3 (b) and (c), the magnitudes of I(t) at both values of V are seen to decrease as the quarantine rate β 2 increases. If the quarantine process of E population is carried out as early as in two days, the disease spreading is greatly controlled (roughly 50 infected population in our simulation of the model in consideration). This highlights the importance of quarantine as early as possible in densely populated places during an epidemic. Hence, the simulation results based on the data indicates that the disease spreading in Delhi is still not controlled properly and needed to take up serious precautions in the state, such as, proper quarantine, social distancing, lockdown etc. . We observed an overall decline in I(t) for all V ≥ 1 [D = N ] leading to proper control of pandemic by monitoring infected population density D as discussed in the above cases. We also found the effect of decreasing the number density n of the population (may be due to policy, such as, social distancing as mentioned earlier) at early times of the epidemic could be a strategy for intervening the pandemic. Again, in Fig. 3 (e) and (f), we studied the impact of quarantine (β 2 ) on the time evolution of I(t) for two different values V = 1 [D = N ] and V = 5 [D = N 5 ]. In Fig. 3 (e) and (f), the magnitudes of I(t) at both values of V decrease as the quarantine rate β 2 increases. In Fig. 3 (e) , if the quarantine process of E population is carried out in ten days, then the I(t) curve becomes flatten, and, if the quarantine process is introduced in two days, then the I(t) curve shows a decreasing trend. Further, in Fig. 3 (f) , if the quarantine process of exposed E population is carried out in two days, then the disease spreading is immensely controlled. This again highlights the importance of quarantine as early as possible during a disease outbreak. Hence, disease spreading in Kerala is quite controlled as compared to other Indian states. . The results showed that for small values of β 2 , the dynamics of I(t) follow Malthusian law, whereas, for significantly large values of β 2 , the curves just start flattening indicating pandemic controlled behavior. As obtained before, if the quarantine process of exposed E population is carried out as early as possible during a disease outbreak, then the number of I(t) can be systematically controlled. But still the condition of Maharastra state is alarming as compared to other state, and proper strategy needed to be taken up. . From the results we observed that if proper quarantine of the exposed population (small values of β 2 0.5) increase in population I(t) is quite large and curve flattening takes long time (100-200 days). But proper quarantine of population E is done (large values of β 2 ≥ 0.5), the dynamics of I(t) became decreased, and pandemic can be controlled. The data based simulation indicates that the scenario of West Bengal is also still alarming as compared to other Indian states.
Disease spreading pattern: We, then, studied the how Covid-19 spread in the parameter space (α, β 2 , I) as shown in Fig. 5 . The dynamics of I(t) for all the five Indian states at one particular value of V = 1000 for thirty realization each show that the trajectories show peaks around 50 days [ Fig. 5 left panel] and start declining. The dynamics show that Covid-19 spreading in Kerala state is quite controlled as compared to other states. Then we calculated I(t) after simulating hundred days in the parameter space (α, β 2 , I) for India and her five states [ Fig. 5 (a)-(f) ]. From the plots, we observed that the population I is relatively large for a large value of α and small value of β 2 . As the value of β 2 increases, the I population drops, and stochastic fluctuations can be seen in these heat maps. Hence, in order to control the Covid-19 pandemic, the parameters α and β 2 need to be optimized.
We have studied the stochastic SEQIR model in the context of Covid-19 disease dynamics in India and its five comparatively worse affected states using stochastic methods. Our numerical simulation results show that policies like social distancing and quarantine have important roles in controlling the disease outbreak and we propose to optimize these two parameters to effectively intervene in the disease transmission. An important consequence of employing a stochastic method is the importance of demographic fluctuations, which are quite visible in the simulation results, to affect the disease dynamics and even in intervening the disease spread. This demographic stochasticity is generally neglected in its deterministic counterpart which is quite important in regulating any system dynamics. Hence, our stochastic simulation method could capture the demographic stochasticity which is non-negligible. We would like to mention that we do not intend to give quantitative predictions here. One limitation of the model under consideration is that by construction, populations from S and E compartments make transitions to the Q compartment where they are assumed to interact homogeneously. This may give rise to a more infected population and we do not see the trends of I(t) converging near zero over relatively less 300 days in our simulation. This also points out that policies such as social distancing and quarantine of the exposed population are not sufficient enough to end the COVID-19 disease outbreak. Other policies like complete lockdown and more testing of susceptible populations should be considered and must be incorporated systematically in mathematical models. 