The sudden outbreak of the pandemic of the Novel Coronavirus (COVID-19) has caught emergency attention for its unprecedented rage in most parts of the globe, overwhelming the healthcare facilities, locking down cities, hitting hard the economy [1, 2] . Within the first four months from the very first report from the city of Wuhan at the end of December 2019 [1, 3] , the epidemic has already reached more than 200 countries, in all the inhabited continents, infecting more nearly 2.5 million people globally [4] . The death toll is already touching the sky at a devastating pace [4] .
Such pandemics, however, are not new in recent times: The 1918 Spanish flu resulted in 50-100 million deaths throughout the globe [5, 6] . In recent years, similar coronavirus outbreak happened in 2003 (SERS-Cov) [7] and 2012 (MERS) [8] killing around a thousand people in China and the Middle east. Following that the outbreak of Ebola cost more than ten thousand lives in Africa [9] . After resulting in massacres in Europe, the pandemic has turned its head towards the United States where more than half of a million people have already been infected within a short span [10] . The complexity in understanding and estimating the spatial spread of the pandemic lies in the conflicting numbers of the rates at which the virus is transmitted to the neighbourhood [11] . Specifically, when multiple hotspots coexist together, only a few recursive events of local transmission turn into a catastrophic avalanche with dramatically increasing number of infections in a very short time window [12] . Moreover, the adaptability of the virus in ever-changing environments is also surprising [13] . Thus, identification of the true character and the self-organization of the epidemic has emerged as one of the biggest global challenges in human history. This is pivotal to develop non-pharmaceutical antidotes with sophisticated intervention strategies to control the dispersal of the pandemic until some clinical solutions arrive in support [14] .
Indeed, there are many major challenges in predicting the spatial spread in a country having multiple epicentres of human transmission, using phenomenologi-cal models [15] . Since an effective model requires the prior knowledge of transmission and recovery rates, it is almost impossible to estimate them in the midway of an epidemic, for the rates depend on the nature of instantaneous social preparedness that vary locally. For instance, Germany had three major infected regions, Bavaria, Baden-Wuerttemberg and North Rhine-Westphalia which had roughly twenty thousand cases in each of these federal states on average by the first week of April [16] . In contrast, the other neighbouring states, Rhineland-Palatinate, Saarland, Hesse had surprisingly smaller numbers [16] . This situation is significantly different from Hubei that alone had more than eighty percent of the total cases in China [17] .
"Essentially, all models are wrong, but some are useful" [18] . Here, we address the dynamics in an epidemic using a model of interacting hotspots. We consider a model state consisting of multiple cities that exchange infected populations pair-wise stochastically where the local dynamics in the individual cities are governed by the deterministic Susceptible-Infected-Recovered (SIR) model with instantaneously or regularly updated populations. Moreover, the order of the random exchange of the infected populations is controlled by a connectivity parameter. We show that tuning the connectivity parameter, it is possible to control the cooperativity of the infection propagation. However, the heterogeneity in the cooperativity determines controls over the outbreak. We validate the qualitative features of the model with the data of the heterogeneous evolution of COVID-19 in different provinces of Germany. We discuss further possibilities of development.
We consider model dynamics in a state consisting of different cities that are allowed to exchange populations among themselves. The population exchange is subject to a containment measure modelled via a connectivity parameter which is assumed to be a constant and identical for all cities. M such interacting cities are considered with initial total populations N i being the sum of the Susceptible (S i ), Infected (I i ) and Recovered (R i ) arXiv:2004.12799v1 [physics.soc-ph] 24 Apr 2020 populations in a city i. For simplicity, the exchange process uses the criteria of classic random exchange of econophysics [19] . During an exchange each of two cities retains an identical fraction of population η and redistributes a random fraction of the infected population pair-wise. The interactions between a pair of cities (i and j) with infected population I i and I j are governed by the exchange rules :
with
The process also ensures the conservation: N i + N j =Ñ i +Ñ j during the exchange processes. However. after t = 0, the model does not use any restrictions to enforce any conservation in the global population. Here, f R is the random fraction (taken as a uniform deviate between 0 and 1) and η is the connectivity parameter. The parameter serves in the same way as the saving propensity does in the context of econophysics models [19] . Now in the city, i, the time dependent population, N i is sum of instantaneous updated susceptible(S i ), infected(Ĩ i ) and recovered (R i ) populations following the SIR equations of motions [15] :
We consider M 2 number of stochastic random exchanges using Eq. (3) (4) between the cities before we update populations using Eqs. (5-7) (See numerical scheme in the Supplementary Information(SI), SI Note.1 ) in each time step. The number of exchanges are optimised in order to ensure significant amount of exchange of infected population among the cities that drives the process. In each of the time steps, we now solve 3M number of coupled differential equations alongside M 2 [20] number of random exchanges to extract the dynamics. We use the time step δt = 1 for the integration throughout the calculation.
We monitor the dynamics of such a state consisting of M = 64, 1024 and 4096 cities with identical initial population (N i = 1 at t = 0) for different ranges of magnitudes of connectivity parameter, η. Here N i , S i , I i , R i represent the fraction of the total population in the state M . Each city has identical initial infected population as I i (0) = 0.0001 [21] . The initial conditions make sure that every city has identical behaviour for η = 1 when the cities are made isolated by some intervention protocol as, for η = 1, the cities retain all the population and does not actively take part in the exchanges. In contrast, they exchange freely for η = 0.
In Fig. 1 , we show the dynamics of the infected (I) and recovered (R) population in the state with M = 64 for η = 1.0, 0.5 and 0.0. For η = 1.0, the cities are isolated completely. Hence, the behavior of the individual cities are identical to the solutions of SIR equations for a case with identical parameters (See SI Note:2 and SI Fig. 1 ). For β = 0.5 and γ = 0.3333, the peak in I appears at t ≈ 50 with a peak height I P = 0.06 while S and R saturates to steady values S S ≈ 0.45 and R S ≈ 0.55 respectively beyond t ≈ 80. This we also see in the upper panels in Fig. 1 for η = 1.0. The peaks in I also appear simultaneously for all the cities t ≈ 50.
For η < 1, as the exchanges of the infected fraction between the cities set in, the behavioural differences in the cities become prominent. In the middle panels of Fig. 1 , we show the evolution of Infected and Recovered population for η = 0.5. We see heterogeneity in the infection propagation: Not only the peaks appear in different times in the different cities, but also the height of the peaks in I i become distinct for each city. The dynamic equilibrium is attained when I i → 0 for all i simultaneously at sufficiently long time after the individual peaks in the infected populations in the cities are observed. These effects become more prominent for lower η. For η = 0, the exchange is completely uncontrolled as individual cities do not have any restriction in retaining the fraction of infected population. Hence, the distributions of both I and R become more diffused as shown in the bottom panels of Fig. 1 . Now for all these cases, we look at the time evolution of the probability distribution of the infected populations, P (I; t) with increasing t. In Fig. 2 , we show the evolution of P (I; t) for different cases of η. As there is a tiny population of infected population in the beginning due to the initial condition, the distribution at t = 0, P (I; 0) ≈ δ(I(0)) for all values of η. For η = 1, as individual cities show a peak at t ≈ 50, the delta peak of P (I) also shifts to higher values of I (Fig. 2.(a.2) ). Since all the cities behave identically for η = 1, there is no spread in the distributions for all t. At higher t when the infected population → 0, the distribution again shifts to smaller values of I( Fig.2(a.3) ).
The situation changes for any η < 1 as the exchanges bring broadening in P (I; t), specifically at t → 50. For η = 0.5, there is a significant broadening in P (I) for 40 < t < 60 [ Fig.2(b) ]. The nature of the tail in P (I) albeit changes with t. We even observe a shift in the mode of the peak in P i [Figs. 2(b.1) and 2(b.2) ]. The length of the tail decreases [ Fig. 2(b. 3)] and disappears at higher t. For η = 0 [ Fig.2(c) ], the broadening in the distributions is more prominent, even at intermediate t[ Fig.2(c.1) ]. At t ≈ 50, the tail in P (I) exists upto I ≈ 0.375. It affirms that there is finite probability to have cities with extremely large proportion of the infected populations of the order of I ≈ 0.375 which is almost 37.5 percent of the total population [ Fig.2(c.2) ]. With increasing t, the tail decreases in I and the peak reaches ≈ 0 as the epidemic mitigates around t ≈ 75.
In order to understand the role of heterogeneity in the infection propagation, we now focus on the behaviour of the moments of P (I) for different values of the connectivity index, η, at different t [ Fig. 3 ]. The mean of P (I), µ = IP (I)dI is shown in Fig. 3(a) . The behaviour of µ doesn't show much heterogeneity, suggesting that the collective peak appears almost at similar times as for individual cities. This doesn't change much as η is varied. However, the height of the collective peak decrease for smaller η around t ≈ 50. The variance of P (I), σ 2 (t) =< I(t) > 2 − < I 2 (t) > show interesting pattern for t ≈ 50 [ Fig. 3(b) ] as η is varied. Quite naturally, σ 2 (t) → 0 for η → 1. As η is decreased, the cooperativity also decreases due to increasingly larger disorder. For small η and η = 0, σ 2 (t) show prominent non-monotonic behaviour with increasing t. However, the non-monotonic dependence vanishes as η → 1.
The shift in µ occurs due to the simultaneous increase in the infected populations while the peak in σ 2 is due to distinct behaviours of the cities. We now look at the behaviour of σ/µ that captures the competition between the two factors. In Fig. 4(a) , we show the phase diagram of σ/µ. σ/µ has prominent non-monotonic dependence near t ≈ 50 for smaller values of η. The peak defines the timescale of heterogeneity in the collective infection 
We now attempt to understand the COVID-19 data of the cases in different provinces of Germany using the insights of the model. The time series data for Germany for all 16 provinces are shown in Supplementary Information (SI Figs. [2] [3] [4] [5] . We now analyse the behaviour of σ/µ computed from the distribution of the infected populations in provinces of Germany (See SI Figs. [2] [3] [4] [5] . We see σ/µ has a sharp peak on 23rd March, 2020. In order to verify the heterogeneity, we further look at the dynamic statistical correlation between the daily new cases in the different provinces for the successive days, using the Pearson correlation between the distribution of the infected populations at all the provinces in time t and t + 1, given by ρ(t, t + 1) = cov[P (I; t), P (I; t + 1)] σ(t)σ(t + 1)
We show the dynamic statistical correlation, ρ(t, t + 1) between the number of cases on successive days from all the federal states of Germany in Fig. 5 . We avoid the double count in the data. The trend in the statistical correlation behaves non-monotonically in the period, with a dip on 23rd March when the peak in σ µ appeared, reiter- ating the effect of heterogeneity due to lack of dynamical statistical correlation. Surprisingly, the date exactly matches the date of the beginning of the lockdown in Germany. From the first week of April, there is a clear trend of flattening of the curve with significantly smaller number of cases with respect to the maximum numbers recorded on 27th March. Always, there are continuous exchanges of populations between the provinces in a country unless any intervention is imposed. The mobility between the cities without restriction, in principle resembles the η = 0 case of the model while the situation with strict intervention behaves as η = 1 case of the model when the cities are considered isolated. Considering time dependence in η, it is possible to realize a dynamic picture from pre to post lockdown situations. The behaviour of σ/µ qualitatively features the balance between the spread of the infection in the cities and the relative change in the collective spread of the epidemic in the state. The extremum in σ/µ captures the extreme heterogeneity in the infection propagation as it optimize the relative balance between the local spread and drive in the collective transmission. On this note, the model that we use delivers the similar qualitative feature of the data of COVID-19 cases in the states of Germany.
Whether the appearance of the peak in heterogeneity is a universal confirmatory signature of containment, can be confirmed only by further studies. There are many possible ways to develop the minimalistic model. It is possible to include time dependence in the parameters. Random deviates of a different kind can also be tested for modelling specific situations. Also, possibilities of further customization needs system specific details and the set of governing equations where specific rules of death, quarantine or hospitalization counts could be considered in addition. Further, it would be challenging to model the transportation of the infected populations using deterministic equations [22] , vis-a-vis the data driven approaches [23, 24] . On this note, we believe, that this study, opens up horizons in exploring the rich dynamics of infection propagation inluding sophisticated intervention among the interacting hotspots.
The author is indebted to D. Syam for continuous encouragements and numerous insightful discussions, P. Das for technical support and P. Chaudhuri for comments on the draft. The data of cases of COVID-19 for the federal states of Germany was taken from the public archive of the situation reports of the Robert Koch Institute on COVID-19 at www.rki.de. There are no conflicts to declare.