Black Death—model and simulation

Accepted Manuscript Title: Black Death − model and simulation Author: Orlando Silva PII: DOI: Reference:

S1877-7503(16)30124-7 http://dx.doi.org/doi:10.1016/j.jocs.2016.08.001 JOCS 534

To appear in: Received date: Revised date: Accepted date:

10-12-2015 30-7-2016 4-8-2016

Please cite this article as: Orlando Silva, Black Death − model and simulation, Journal of Computational Science http://dx.doi.org/10.1016/j.jocs.2016.08.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Black Death – model and simulation Orlando Silva Universidade de Évora, Portugal

Highlights:    

One continuous deterministic SI epidemic model is chosen and further developed to simulate the Black Death. A model for the interaction between short range dislocations and contagion is proposed that minimizes the need of a classical diffusion term for infected transport. One modified expression for the classical contagion rate law takes in account the saturation effect in contact rate. The results of the spatiotemporal simulation of Black Death are compared with historical data.

ABSTRACT The Black Death analysis and data organizing have been performed in unifying manners, among others, by Byrne and Christakos. The mathematical modelling of Noble of 1974 did not originate full spatiotemporal 2D simulations. The present model incorporates a new term to calculate the effects of the interaction between contagion and short range dislocations. The contagion rate term is modified to account for saturation effects in the contact rate. The model and parameters used very nearly reproduce important features of the spatiotemporal evolution of the Black Death, the calculated global mortality attaining a very large fraction of the European population. Keywords: Black Death, epidemiology, model, simulation

1. Introduction Approximately at the autumn of 1347, one disease was introduced in Europe via sea ports, possibly by land too, that spread to all the continent by sea coast, river and land routes, resulting in the infection of almost all European regions. It came to extinction by 1351, meanwhile killing a large part of the inhabitants, usually estimated between one fourth and half the population. Attributed to the Italian notary Gabrielle de Mussis is the claim that the disease was transported by boat after the infection of Caffa, in Crimea, till Messina and then to Italy. Differently, Wheelis [31] discards the role of Caffa in the contagion of the European ports. In the following decades, some other epidemics led to much lower mortalities in Europe, possibly by the rise of the fraction of the immune population or due to the change of the pathogenic agent. Most authors support the disease was caused by the bacteria Yersinia pestis while others reject this hypothesis. More detailed information can be found in [13] and references there. According to [4] and most authors, the transfer of the pathogen to humans is made by infected flees from rats. Others sustain the transfer is made by infected fleas and lice between humans (see [16]), the goods in trade routes possibly carrying infected fleas even in the absence of rats.

1

Other modes of contagion has been mentioned, as the infrequent transfer of pathogenic agent from organic material from soil, [18,21]. It is possible that all these and other transfer modes coexist in different proportions, some of the modes possibly changing characteristics in space and time. A great number of historical reports and studies of the Black Death exist. Qualitative descriptions and quantitative analysis (for example [10,13]) very much contribute to more precise knowledge about the facts and the evolution and effects of this epidemic. Some studies of the spatiotemporal dynamics of the Black Death lead to quantitative simulations of the epidemic spreading, including the effects in European population, Noble [23]. One 2D spatiotemporal simulation of the Black Death was however not performed in the study of Noble. In the present study one mathematical model of the same type of Noble’s one is built. The phenomenon of interaction between contagion and pendular population movements is added to the model. The calculation of the contagion rate is made by means of a modified mass action term that incorporates the effect of saturation in contact rates between humans. The performed calculations give rise to simulations of the spatiotemporal evolution of the Black Death. The precise way in which an epidemic propagates, along with its effects in populations, may involve a very large quantity of elements and interacting events, the considered most important ones being included in the model. Pathogen transfer between infected and susceptible individuals, births, deaths and populations’ motions, are some of the important phenomena to be included in an accurate description of the epidemics evolution. The dynamical model of the Black Death can then calculate the total variation with time of infected and susceptible at each location, in terms of the entry and exit fluxes of individuals, of the rates of change from susceptible to infected, and of mortality rates. For the Black Death, the direct historical data provides limited information about the rates and fluxes at each place and time. The spatiotemporal evolution of the populations under epidemic situation is thus not accurately described for the whole space and time domains. Better descriptions of the epidemic evolution so need, for each epidemic phenomena, models that quantify the insufficiently known terms, then including them in balance equations of each one of the epidemic species. The details about the contagion modes and the involved biological species must be known, or in some way assumed, in order to calculate the rates of infection of the biological species of interest. If the cross contagion between different biological species is frequent, there will be a close relation between the wave fronts of infected of the species transporting the agent. If cross contagion is very rare, the wave fronts of the infected of different biological species will be almost independent. It is noticed that the spreading of the disease only by rats would occur at much lower velocities than the spreading velocities in the case of the Black Death, [16]. If human infection propagation depends chiefly on rats’ unperturbed infection wave, the slow infective wave in rat populations could determine low velocity of human epidemics, unless rat infective wave can be made quicker traveling with humans, what would make both wave fronts to approximately coincide. Several contagion modes can coexist and can in principle be modelled. Under the uncertainty about the more important infection modes, some important simplifications will be made, leading to a model only dependent on the balance equations for susceptible and infected populations of the human species. The model is supplied with values of the coefficients and with initial conditions that are supposed to approximately describe the initial epidemic facts and properties.

2

The resulting differential deterministic equations will be transformed in algebraic equations, whose solution will be obtained by numerical methods. The numerical results will be compared with available data of mainly historical origin, after the work of Christakos et al [13], which in several occasions will be designated in a simplified manner by [13]data. These very valuable data are built based on historical qualitative and quantitative data, on hypothesis about some characteristics of the infective agent, the agents of transport and transfer, population distribution, geographical and climatic environments, and a multiplicity of other pieces of knowledge, scientifically established or facts related to the particular epidemic. From all these sources, stochastic methods were used to generate more consistent information from fragmentary data. The dates of beginning and end of epidemic at one place, rates of contagion, population size, mortality including calculated time distribution mortality based on global data, are some of the data directly acquired or determined by methods there described. The very important effort dedicated to gather multidisciplinary information, to the development of methods and to their application results in a coherent and well-structured set of results that will be here used as basic data for our work, and data for comparison with our own results. Due to that dependence, the quality of our results is generally expected not to overtake the quality of their results. Our calculated results, sometimes termed C-results, will tend to smooth irregularities present in [13]-data. In some situations, more research will be needed to decide if regularity or variability are the right picture of the facts. The [13]-data about the beginning of the pest in [13] do not explicitly display one uniform criterion to the pest start. In many cases they will surely be the dates of first deaths or when symptoms are clear and death follows. In a fraction of cases, especially when news about pest are already abundant, the arrival of voyageurs with information could be taken as pest beginning, be it the case or not. The uncertainty about the criteria of pest start in many of the [13]-data raises the problem of defining how to compare with C-results and what criteria to use for the beginning of infection in entrance ports and inland infection. These quantities can further assume different values in the European spatial domain, probably also in the time domain. Dealing with the weak reliability of contemporary sources, [13] cross-validate and filter the dispersed and disparate information as to categorize and treat it to extract the most acceptable data characteristics, quantified when possible. Even so, the referred dates are not enough unambiguous. Some practical options about criteria about the dates of local pest start, not sufficiently solid however, will be presented later.

2. Model The evolution of one epidemic system depends on a large number of elements and on the relations between them. In the case of the Black Death epidemic important knowledge limitations exist however. The multidisciplinary study [13] brings together important disperse information to which stochastic methods are applied to reconstitute some of the important missing information. Quantification of miscellaneous qualitative and quantitative information extraordinarily improve the accessibility and uniformity of data. Regularities and variability that were not easily found in the dispersed original information about Black Death are now easier to seek. The model to be applied in the present work is intended to reproduce some of the main characteristics of the Black Death, avoiding many of the system finer details. The network of long range travel routes with their ramifications practically extends to all European regions carrying goods and people, eventually carrying the infective agent. A large number of short range, fixed or varying, dislocation paths further increases the capability to reach all the European territory. The accessible European spatial domain will be considered as a continuous two dimensional space, where all the fields of variables and coefficients are defined. The low fraction of 3

urbanization (with reference to cities with more than 5000 inhabitants), around 10% as indicated by [5], further supports the continuous model option. The systematized data suggests a significant regularity in the spread of the epidemic initiated at seaports (Fig. 3). The observed irregularities of the spreading velocity are believed to be related either to natural and political barriers, as to existing lines of rapid progression like long range trade routes. Some variations of the velocity can be related to weather conditions dependent on the seasons of the year. The fact that the overall spreading hits almost all European regions, at considerable velocities, reinforces the belief of significant efficiency in transportation and contagion of the pathogen. These regularities allow a deterministic description of the epidemic system, originating a continuous and deterministic mathematical model of the epidemic system. A SI epidemic model (see for example [22,23]) is considered to be appropriate to the description of this epidemic system. The mortality of the infected is very high and the possible survivors are thought to have small influence on the epidemic evolution. The rates of change of the densities S and I at each location and time can be expressed in terms of the rates of contagion, births and deaths, the fluxes of S and I resulting from long range dislocations, and the rates of change of S and I as result of the interaction between short range movements and contagion. Adding these rates and transport fluxes results in the total rates of variation of infected and susceptible populations at each location and time. The equations so obtained represent the dynamics of the epidemic system as expressed by the mathematical form of the model. The initial values of variables and parameters are prescribed on the basis of the available data, adjusted for better fit between data mainly of historical origin and the simulation results using the present model. At times and places at which infections are believed to come from sea, infective sources are imposed, in this way originating epidemic waves. The balance equations of the model in differential form will then be approximated by algebraic equations, obtained by the finite volume method, see details in [28] and origins of the method in [3]. The solution of these equations leads to the spatiotemporal evolution of the system as described by the model. The results will be compared with data of multidisciplinary sources, particularly historical, analysed by methods of [13,14], which here receive the label of [13]-data. The comparison between numerical results and treated historical data are expected to have the capability of suggesting, besides the already performed, further alterations to model form, values of the parameters and particular conditions, and possibly to question some of the historical data. Next subsections include the submodels that allow to calculate the net rates and fluxes that contribute to the rates of change of local populations defined in terms of their densities S and I. 2.1. Creation and destruction rates of susceptible and infected Creation and destruction rates of susceptible and infected by effect of the contagion phenomena, expressed in terms of their densities, are modelled, [22], by means of 𝜕𝐼 | 𝜕𝑡 𝑐𝑜𝑛𝑡𝑎𝑔𝑖𝑜𝑛

= (𝑘𝑐 𝑆)𝐼

(1)

𝜕𝑆 | 𝜕𝑡 𝑐𝑜𝑛𝑡𝑎𝑔𝑖𝑜𝑛

= −(𝑘𝑐 𝑆)𝐼

(2)

The contagion is active after the latent period, being however considered that all infected are also infectious. 4

This model for the contagion rate has its remote origin in chemical reaction studies, in which kc normally varies with temperature, pressure and the presence of third bodies, so inducing reaction rates to vary. The applications of this model extend to biology and ecology. Independently of its correctness as applied to chemical processes, its application to epidemics must be performed with extra care, making attention to the influence of several conditions that may cause variations in the value of kc. One such influence can be conceived in respect to the effect of the density of susceptible in equations 1 and 2. The infected person, like any other person, will not go into two bakeries in its way home, and will normally not visit all farms when coming from the farm where has already worked for several hours. In this way, the increase of S in one particular zone does not lead to one proportional number of interactions per day of one infected with all the susceptible with who the infected could contact. Other voluntary or involuntary contacts may not be subjected to one saturation phenomenon, the contagion rate keeping the proportionality to S. One detailed model for the contagion parameter may involve an arbitrary amount of details and will not be addressed here. To include the effect of the saturation of the contact rate, one ad hoc expression is devised to the contagion parameter: 𝑘𝑐 = 𝑘𝑐0

𝑆 𝑘2 +𝛼𝑆 2

(3)

The value k2=1 km-2 has been chosen. The limit of the product kc.S is kc0/ as S→∞. Using the values of Table 1, the rate of contagion is roughly proportional to S at values S < 15km-2, for larger S values approaching one saturation rate only slightly dependent on S. The proportionality with I is maintained as shown by equations 1-2. [13] state the duration of epidemic in cities to be proportional to the number of inhabitants. This is

compatible with equation 1 if the city is not considered one well-stirred reactor but one spatially extended system where the epidemic wave spreads at limited velocity. In this respect the expression of equation 3, or another of similar type, may be important in limiting the velocity of spread. The complete equations 1 and 2 would need to be solved inside cities, what will not be done here. The death rate of the infected as consequence of disease may be approximately modelled by 𝜕𝐼 | 𝜕𝑡 𝑑𝑒𝑎𝑡ℎ

= −µ𝐼

(4)

A death rate term of the susceptible population is included, during a time interval tc (see Table 1). This term calculates the death rate of susceptible due to augmented mortality in the situation of economic and social crisis after the epidemic wave arrival: 𝜕𝑆 | 𝜕𝑡 𝑐𝑟𝑖𝑠𝑖𝑠

= −µ𝑐𝑟 𝑆

(5)

One similar term for the infected will not be included, as death rate by disease is assumed to be much higher than this type of contribution. Also the birth rates are not included because they have a minor contribution to the total variation rates for this particular epidemic and normally are close to equilibrium with the natural death rates. 2.2. Transport fluxes of susceptible and infected at long distances The pathogen transport by the infected, particularly at trade routes, is essential in the epidemic spreading and the fluxes must be modelled and calculated.

5

Caravan and boat transportations can be viewed, in average terms, as a moving continuous media where infected, pathogen and transport device move together. In this way, the infected fluxes can in principle be calculated in a convective form. Nevertheless, details about this transport, like the frequency of the events, the dates, directions and distances reached by these transports, densities of the infected in the streams, are far from known. An alternative way to model these transport fluxes is so highly recommended. A gradient dependent diffusion can approximately reproduce the transport of infected, [22]: 𝜕𝐼 | 𝜕𝑡 𝑟𝑜𝑢𝑡𝑒

= ∇. 𝐷∇𝐼

(6)

The diffusion coefficient D can be estimated from convection data, if available, by means of 𝐷=

𝑓𝑑 2 2𝜏

(7)

where f is the fraction of the population at the departure location that is included in the crew, d is the average distance reached by the expeditions that take place in average every  time units. The approximate knowledge of population consumption rate of outside goods, the carrying capacity of transport devices, the average travel distance of the trade excursions, the population fraction involved in long distance trade, the frequency of travel expeditions, would allow estimate D coefficient. In the present study, a rough estimation of this parameter was made, in which, besides consumption, other travel events may also be present, like military, political and religious expeditions. The resulting D value was replaced by a reasonably close value (Table 1), leading to a better adjustment between simulation results and historical data. In the diffusive formulation, the transport of infected will always be oriented in the opposing direction to the gradient of infected density. The infected flux oriented to the lower infected density has the effect of spreading the epidemic in the way normally happening in reality. Moreover, in Black Death, the lifetime of infected is presumably shorter than the average time of long range two way travel. So there is no need to calculate a counter gradient return flux, suggesting the diffusive formulation to be appropriate. In cases in which the infected move to places of higher density of infected, as can occur in reality, the infections are already taking place in the destiny place, so not significantly modifying the course of epidemic reality. So, the effects of pro-gradient infected flux are usually not very important. In fact, the value of the infection rate exceeds the diffusion flux after the time 𝑡=

1 𝐷

𝑘𝑐 𝑆+(𝛿𝑥)2

(8)

where S is the susceptible density in the destination place and x the distance between the place of flux origin and the place to be infected. For the conditions of Black Death, the estimated value of the time t at which the diffusion flux ceases to dominate infection rate is very short. Then, once one place is infected, very soon the local evolution of the system is dominated by the source terms and no more by the value of the transport flux. In the present model, the diffusion coefficient will have a small constant value at the majority of European domain, and large values at trade routes that occupy a very small fraction of the European area (Fig. 1a). It will exhibit very small values at natural and some political or cultural barriers (Fig. 1b). In the present model, the flux of susceptible is supposed to have a minor demographic effect if compared with the rate of transformation of susceptible in infected. For this reason, no term for the 6

diffusive transport of susceptible is included. The effect of introducing a diffusion term in S equation would lead to the, supposed historically not too important, smearing of susceptible by the effect of their transport from the high to the low S locations. 2.3. Interaction between contagion and short range motion The infectious agent transportation outside the main trade routes is mostly related to the short range population movements in short range trade routes and in the dislocations of the everyday life. In their daily occupations, a fraction of the population displaces from their habitat to nearby places and come back after not long time. These pendular displacements are executed by infected, that may infect susceptible at the neighbourhoods, and by susceptible, that may return in the new condition of infected. The interaction of pendular motion with contagion is modelled by estimating the fraction of dislocated population and calculating the infection rates resulting from these types of dislocations. During a fraction of time fIt, a fraction fIp of the local infected will leave their place, so leading to the absence of local contagion caused by them. A number of infected will come from neighbourhoods and will infect the residents at a certain rate. The fraction of dislocated varies in an unknown way with the distance r. A convenient potential function is used to represent the fraction of infected that travel till certain distances, 𝑎 𝑏

𝑓𝑝𝐼 = ( 𝑟 )

(9)

The exponent will receive the value b=2, while a=0 is taken for large r values. Naming Ipend=ab.fIt, the net local infection rate caused by these movements may then be modelled at small distances by 𝜕𝐼 𝐼 | 𝜕𝑡 𝑝𝑒𝑛𝑑

𝐼 = 𝑘𝑐 𝑆∇. 𝜎𝑝𝑒𝑛𝑑 ∇𝐼

(10)

where the upper I is used to identify the rate of infection due to the motion of infected. The pendular motion of the fraction fSp of susceptible, during a certain time fraction fSt, implies that they will not be locally infected during that time, but can be infected in the neighbourhoods, at a rate that depends on the neighbourhood infected density. Some of them return home in their new condition of infected. The net local infection rate originated by the susceptible pendular motions may be modelled by 𝜕𝐼 𝑆 | 𝜕𝑡 𝑝𝑒𝑛𝑑

𝑆 = 𝑆∇. 𝜎𝑝𝑒𝑛𝑑 ∇𝑘𝑐 𝐼

(11)

In this epidemic, at which the infected live during about a month and have symptoms during less than one week, we will consider Ipend≈Spend , accounting for the relation between the time fraction of displacement of the infected, close to the time fraction for the susceptible. If Black Death would be a bubonic plague, Ipend≈1/2 Spend would be estimated, to account for the relation with the time fraction of absence of the infected, expected to be approximately half the time fraction for the susceptible. 2.4. Model equations of the epidemic dynamics Balance equations for infected and susceptible accounts for the source rates and transport fluxes that determine the variation rates of each epidemic species S and I of one biological species.

7

Simplifying the terms of interaction between contagion and pendular motion, the balance equations for I and S densities are 𝜕𝐼 𝜕𝑡

= 𝑘𝑐 𝑆𝐼 − 𝜇𝐼 + ∇. 𝐷∇𝐼 + 𝐷𝑝𝑒𝑛𝑑 ∇2 𝐼

(12)

𝜕𝑆 𝜕𝑡

= −𝑘𝑐 𝑆𝐼 − 𝜇𝑐𝑟 𝑆 − 𝐷𝑝𝑒𝑛𝑑 ∇2 𝐼

(13)

𝐷𝑝𝑒𝑛𝑑 = 𝑘𝑐 𝑆𝜎𝑝𝑒𝑛𝑑

(14)

pend is the sum of the homologous parameters defined above for susceptible and infected. Analytic solutions of equations 12 and 13 do not exist, the finite volume method being used to numerically calculate the spatiotemporal evolution of the epidemic system subject to the rules of the chosen model. The algebraic equations are solved in time by one implicit formulation or by the 4th order Runge-Kutta method, the results being very similar when using a time step of 0.2 day. 2.5. Biologic species contributing to the epidemic and effects in model building Model equations 12 and 13 describe the dynamics of the populations of one unique biological species. More non-parasitic species may however be active in the epidemic dynamics. Different hypothesis exist about the number of intermediate vector species and the modes by which the infective agent is carried and reach humans, see for example [19,9]. Several parallel modes can even be simultaneously active with different intensities. Several years after the description by [26], the international scientific community came to accept the infection of humans to be caused by infected fleas transferred from rats [10]. When large mortality of rats occurs, the infected fleas transfer to humans, and the bacteria transferred to humans then spread among them. Scott and Duncan [25] however sustain the contagion between humans to be practically non-existent. The system dynamics modelled only by equations 12 and 13 can be extended to include rodent populations and even fleas populations, [19]. The model equations must then include one or more cross-contagion terms to calculate the effects of human contagion by other biological species. If the rates of cross contagion are similar or larger than the intra-species contagion rates, the wave fronts for rats and humans are expected to approximately coincide. The contagion rate can depend on similar form on the epidemic species as in equations 12, 13 and 3, but the contagion coefficients may need modifications, including variations in space and time, to accommodate rat contact rates and the effect of cross contagion with a rat population of unknown density. The two equation model can so be maintained at the price of imprecise kc. In Simond contagion mode, the average contagion rate from rats to humans is very small. The spreading of the rats’ infective wave is very slow as referred by [16], unless humans transport rats and fleas in the trade dislocations luggage. When staying at locations with high rats’ density, according to Simond mode, fleas abandon humans. The effects of this transfer in the epidemic spreading must then be investigated. The regular epidemic spreading could halt for long periods in some locations, waiting for high rat mortality at those places. The inclusion of balance equations for rats’ populations would be needed. Historical data show no evident pauses at point locations, although they do happen when crossing certain lines, apparently related to barriers to human dislocations. The lack of evidence of large rat mortality during the Black Death, [12], suggests this contagion mode not to be the most important one during the epidemic spread, excepting in the supposed rare cases in which great mortality of rats happen. 8

Other authors, like [16], state that the contagion to humans is essentially done by the transfer of infected fleas and lice between humans. [12] refers that no epizootic of rodents has been reported immediately preceding or accompanying the human plague. If this is the case, a model with a unique biological species is adequate to the dynamics of the Black Death, and the balance equations 12, 13 for human epidemic species will suffice. If the disease spreading is fast as in the case of the Black Death, it is probable that this contagion mode exists, as propagation velocity by rats alone is slow. [16] and several references there point to the adoption of the human species as the main carrier of the infective agent. [14] and references there support the idea that Black Death must not be a bubonic plague but any other type of epidemic. This one will be the adopted version in the present work, with the effects it may have on the model and respective parameters, particularly on mortality coefficient. Other modes can be active in parallel, also the pathogenic agent may not be unique, and this case may be not definitively closed. The simplest contagion mode is included in the present model, compatible with Hufthammer mode. The human contagion is considered to occur between infected and susceptible humans, its average rate being supposed independent of the numbers of fleas and lice. 2.6. Model specifications In order to use the mathematical model of section 2.4 in simulations of the Black Death, some quantitative specifications must be made. The values of the parameters that enter the equations; special conditions like initial population distribution and the locations and times at which special infection sources are activated that correspond to infections from the outside of the spatial domain; density values for the infected at these harbours; threshold values for I that are sufficient for the epidemic to maintain; and criteria for epidemic perception – all these and possible others look for quantification that can lead simulations to supposed realistic results. Values of main parameters

Historical, medical, biologic, climatic and geographical data help in the specification of the parameter fields. In the absence of detailed knowledge about several coefficients, model specifications were built on the basis of initial semi-empirical values. These values are used in simulations and the results compared with [13]-data in order to choose parameter values, Table 1, that optimize numerical results.

The basic coefficient values used in the model equations 12 and 13 are shown in Table 1. Routes, barriers and time dependent transport flux

The observation of the dates of infection of cities from [13]-data suggests the existence of large or small epidemic wave velocities, according to the seasons of the year, possibly also to the existing travel routes and to a small number of barriers that alter the regional velocity of plague progression. Lines of high and of low diffusion coefficients, where initial population density also have similar modifications, are then defined in the otherwise uniform diffusion fields. Geographic, political and cultural barriers may delay the epidemic progression. This effect is introduced in the model parameters by lowering, with different factors, the diffusion coefficients for the most important barriers. The low diffusion lines, depicted in Fig. 2, result in delays of the calculated epidemic waves. A few of those lines were not related to known barriers, but imposed on the basis of the delays existing in the progress of the patterns revealed by the [13]-data.

9

Some lines with high diffusion coefficients are defined to reflect and simplify the high velocity paths associated to the system of trade routes. The system of high diffusion coefficient, to be used in the simulations of the present model, must be considered a crude simplification. As one example, the only one high diffusion line in Great Britain (Fig. 1a) must be considered as one synthesis of the long range trade route system of the island (see [15] and references there). The lines of high diffusion coefficient partially coincide with the system of long European valleys and rivers. Some high diffusion lines parallel to the coast are supposed to replace the effects of coastal shipping. The velocity of the epidemic waves is a raising function of the diffusion and contagion coefficients and of the initial susceptible density, but a decreasing function of  (see for example [22] for the case of a similar mathematical model). The month mortality also depends on these coefficients. When searching for good fit between C-results and [13]-results, D, kC, S0 and , may vary within moderate intervals, still leading to similar results of the propagation velocity.

Trade, agricultural work and general dislocations are reduced in winter, so the need to introduce an ad hoc time transformation to the diffusion coefficients, that takes the form 𝐷 = 𝐷0 . (1 + 𝐴. |𝑠𝑖𝑛 ∝|. 𝑠𝑖𝑛 ∝)

;

2𝜋

∝= 365 ∙ (𝑡 + 165)

(15)

where t is the time in days after the beginning of October 1347, and A=0.33. The same transformation is applied to the values of Dpend in equations 12 and 13. Different specifications could be appropriate at different regions, maybe at Scandinavia, but will not be applied here. The value 165 of the phase implies the maxima diffusion coefficients to occur around June, what coincides with the conclusions of [8]. It must be noted that the ratio between maximum and minimum values of the diffusion coefficients is not directly reflected in Börner coefficients. In their calculations of the velocity of the epidemic wave front, a number of other terms are included which, in our model, must be considered as implicit in the diffusion coefficients, making the two forms of calculation not directly comparable. In a number of cases, aggregates of cities displaying simple describable geometries are associated with apparently sudden delay of the epidemic spreading that affect relatively large areas or considerable number of cities. These events are normally interpreted as consequence of coherent effects of some type of barriers to transport so affecting diffusion coefficient and population density along the defined lines. The barriers are mostly related to geographic obstacles, possibly also political, administrative and religious frontiers. Some delays may not be independent of rigorous winters. The lines of apparently quicker spreading normally coincide with rivers or routes active at that time. The definition of diffusion coefficients dependent of time is related to the different velocities of Black Death at the year seasons. The case of Scandinavia would maybe require a different time dependence of the diffusion parameter. Another tool to adjust regional epidemic velocity is related to the definition of the initial population density, although commonly respecting the existing population data, namely [2]. The irregularities observed at relatively small space scales are commonly not subject to any action on parameter values.

10

The definitions of minimum densities of infected, I for perception of infection and initial I at harbours infected from sea, were determined by mixing guess, empiricism, simple scale analysis, all of them suffering limited optimization to better fit of C-results to [13]-data. Initial population density field

The initial susceptible density field was built by freely adapting a number of published data, not all of them referring to the same epoch, [2,11,17,24]. These freely integrated data were submitted to later corrections as effect of the comparison of C-results with [13]-data. Fig. 2 is so a coarse map of the population density at the beginning of the Black Death. The densities at vast eastern and south-eastern regions have weaker reliability, due to lack of knowledge of demographic and epidemic historical data. At the rest of Europe, also the used initial S field is not in principle very accurate. In fact, some other values of the ensemble of parameters, initial conditions and particular options may in principle lead to similar results.

Dates of infections and initial number of infected in harbours

In their arrival to harbours, the ships that are carrying the pathogenic agent disembark people and goods, and will be supposed to transfer an average of 3 infected individuals, possibly one too modest figure. These previously infected individuals are the sources of the spreading infective waves propagating inland with origin at some European harbours at specified dates. The dates used in the simulations for these harbour infections very closely follow the ones in [13]. The infection of sea ports originated by the contagion associated to long distance boat travels can be identified in some cases in Fig. 3. In the case of ports infected from sea, the date of the perception of the presence of the epidemic will be normally assumed as the date directly coming from [13]-data. Threshold density for infection spread

If in October 1347 just one infected arrived at one uninhabited large region of Europe, maybe at north Scandinavia, the pest would not spread after this point. But if the infected landed at a small village, pest spread would depend on the communications of the village with other places. The number of diverse situations that lead to the average propagation or extinction of the pest is large. This suggests the need for one separate study to decide the minimum density that in average causes the epidemic to spread in the scope of one deterministic continuous model. In an ad hoc decision, this number is chosen to be I=10-5km-2, roughly equivalent to one infected in (300km)2. Below that number, the density is turned to zero. Criteria for pest start perception

In the ports, the infections coming from the sea will be normally soon identified. In addition to the immediate signs of disease that can exist, the stories of the sailors will persuade the locals that anything, identified as pest or not clearly identified, is occurring. This perception may have one large probability of being transmitted to written documents, locally or destined to inform someone in one neighbour city. So the time of perception will in these cases stay close to the real infection time. In relation with some [13]-data, epidemic start is directly identified by the clear rise of the number of deaths in one city, or indirectly from dates of some other events as shown in [13]. The infection of one city could be coincident with the death of one outsider that at the end of disease period arrived to the city at that date. But the infected outsider could arrive at the beginning 11

of the incubation period, affecting the time of epidemic perception by at most one month. In certain cases however, the incoming infected may stay and infect others. If the city was not aware of the pest, something strange in the overall situation could be perceived maybe only after the second death. How this information would come to be recorded is also uncertain. Other possibilities exist, for example that infected but not ill foreigners stayed for only a short time and then leaved. In average there will be some delay in pest perception in relation to the real infection start. This delay affects local people, chroniclers, historians, and may to a certain extent remain in [13]-data obtained by models applied to original data from miscellaneous sources. Admitting that the criterion of pest beginning in [13] is influenced by the reports of the epoch, and that the information in reports depend on facts like the referred above about the number of deaths, the infective density can be a good indicator. To the local perception, the generally unknown ways in which facts influence reports, rise the uncertainty. As the miscellaneous sources are very dissimilar in many aspects, the criterion to choose the date of pest start is largely arbitrary. In many cases in which the information is scarce and pest arrival is not very clear, the criterion of pest perception time defines the moment at which the number of local previously unexpected deaths is almost two in a small region of about (30km)2. This will correspond to one infected density of about 2x10-3 km-2, which will be the value at which simulations make time comparisons of Cresults with data directly taken from [13], with exceptions in one very small number of cases. This author however is not devoted to the belief that this reasoning is the most adequate to be used to choose one common criterion for the pest declaration in [13]-data and in simulations. In pest importing harbours, there is a close correspondence between the times of infection for both inhabitants and model user. Simulated pest beginnings then occur at the time identified in [13]-data, supposing chronicler was well informed. As stated before, in pest importing harbours, at these times, 3 infected individuals are added. 2.7. Possible improvements to the model Two connected space domains with different diffusion coefficients

It is expectable that the flux of people, goods and infective agents will have different values at different trade routes, between these and countryside and also at different locations in the countryside. It is possible to define one normal 2D domain for the epidemic dynamics at countryside and one 1D domain for the net of trade routes. Between these two sub-domains the diffusive coefficients may be made small so to guarantee the expected effect of quicker infection in the 1D domain with no excessive loss of infected to the 2D domain of the countryside. This would have the expected effect in the economy of the epidemic system. One further 2D domain can be inserted in the interrelated subdomains in one analogous way to reproduce the epidemic dynamics inside large cities. In these new subdomains, the conclusion of [13] about the approximated proportionality between pest duration and city population could be checked. Even the extraordinary possibility of using the same equations and even possibly the same parameters in one more complete model could be tested. Givry [13]-data [13] and references there about daily and monthly mortality disagrees with the general characteristics of monthly C-mortality at general locations. Using the mean-field formulation of present model, monthly mortality results wider in time and more intense, leading to a larger population depletion. The study of the epidemic dynamics inside cities could possibly lead to alterations of the calculate time evolution in cities, different from what is expected using our model defined by equations 12-13 informed with parameter values of Table 1.

12

The relation between the models and the parameters of the two realities in very different spatial dimensions would be useful in deciding the limitations of models or alternatively to check the possible uniformities and heterogeneities of used parameters. The effects of downsizing in the measure of data is referred in [13] and references there, nevertheless mentioning the lack of data at the smaller scales. Effect of variability in non-linear terms

Problems like the effects of time and space variability of quantities inserted in non-linear terms can be partially solved by means of the expression of the former mean-field variables and parameters as the sum of their respective average values and fluctuations. The inclusion of these terms in the equations and averaging then leads to covariance terms of unknown magnitude, so requiring to be modelled because of the non-closed condition of the newly generated equations. When the unknown new terms are modelled in forms dependent of known variables, then the new equations are closed and can normally be solved. Equations of this type normally include not only the rates and fluxes of the previous equations, but also terms that insert the effects of non-linear variations in the general balance equations. Examples of these procedures are found in a different field of knowledge, for example in [30]. The modelling of some unknown terms may be executed with the help of stochastic methods like in [30,13].

3. Results The Black Death simulations are performed using the model here presented, supplied with the parameter values of Table 1 and equation 15, the initial susceptible density field of Fig. 2, and the special conditions imposed to the infected density. In particular, non-zero densities are imposed at specified positions and dates. 3.1 Epidemic spreading In the last months of 1347, Sicily seaports, Genoa, Marseille are reported to have received infection by sea. The spreading of Black Death continued by land, rivers and other seaports by shore and by long distance vessels transporting the pathogenic agent. The epidemic waves progress along seashores as well as from the coast locations towards inland. Pictures of the spatiotemporal simulation results and from [13] data are shown in Fig. 3. The parameter m identifies the months since the pest beginning at October 1347. Epidemic spread in the areas of nowadays countries is outlined next in a loose time order. Some topics about the progress of the historical and calculated wave fronts are mentioned, although in a way far from systematic and extensive. The referred historical data is almost exclusively based on the work of [13]. These data are here commonly labelled historical or [13]-data, although they were determined from multidisciplinary sources using methods described in [13]. Italy

Sicily was infected in October 1347 while the Italic peninsula, according to [13]-data, received the pest almost simultaneously at south and north, close to the end of the same year. Although Genova was infected in December 1347, no quick spread is mentioned to occur to neighbour locations. A spreading front having its origin in Pisa after the initial months of 1348 propagates till Veneto and Lazio regions. Possibly the wave front that reached Bologna propagated then to NW and arrives Torino in November of 1348, the last large Italian city known to be infected.

13

According to [10], the plague would need five to eight weeks to establish in a city, so the reports of appearance may have to be adjusted backward to ascertain the date of pathogen’s arrival. The procedure would naturally apply to large and small cities, possibly with different corrections. In the present study, section 2.6, instead of backward adjustment of [13]-dates, one delay in the declaration of infected place was made for I > Iperception, one I-value thought to correspond to the infected density at cities in which the epidemic was already recognized as having started. Historical data led [13] to choose December 1347 as the infection date of Venice. The surrounding places were nevertheless hit many months later by an infective wave apparently coming from Pisa. A very effective isolation of Venice would be compatible with this description. In this work, a risky exercise was executed, by arbitrarily choosing February 1348 as the infection month of Venice. This was conjugated with the insertion of a barrier of low permeability to transport around Venice, so that good fit with historical data could be achieved. Historical analysis must decide the geographical origin of the infection of places around Venice and the possible lack of easy communications with nearby localities. The comparison of numerical results with [13]-data seems to reveal an early attribution of infection dates to Venice and other large cities. Possibly small localities should have larger imposed delays in the historic dates of infection than large cities. This option could correspond to initial infections with smaller number of infected than the number inserted in large cities. For cities infected from sea, large or small, it is very probable the insertion of one significant number of infected at disembark, resulting in the immediate perception of pest beginning. [13]-data indicates the possibility that all the mentioned Italian cities were already infected before 1349. The coefficients and initial S values used in simulations lead to similar results, the whole Italy being infected at dates close to [13]-ones. For southern Italian regions, there is scarce data and, in the centre-eastern regions, the data of [13] do not mention any infected place. The accuracy of the numerical results cannot then be adequately tested there. The dates of infection of a large number of European cities is not known, although their infection is assured. In cases like these, the evolution of the calculated infective waves allow estimations of those dates. Only the year of 1348 is known as the infection date of Milan, calculated results suggesting a date close to September 1348. However, only in cases of a large number of cities nearby with identified infection dates would this type of estimates be considered reliable.

Balkans, Greece, Turkey

A strong probability exists that these regions have been infected before October 1347, [27]. The historical data suggest the infection of Adriatic east coast by sea or seashore. The scarcity of data for these regions allows a high degree of liberty in the choice of initial conditions for population density and infected harbours. The confidence on the simulation results is therefore quite low. France

After the infection of Marseille in December 1347, the wave fronts evolve to north, to west, to east by seashore, and southwest towards Spain. The western front apparently do not cross the Pyrenees, the numerical results matching this hypothesis. By June 1348 the north coast is infected by sea, the infective wave mostly moving in 14

southern directions. Northward and southward waves join about November of the same year. It is not clear from [13]-data and C-results if Low-Normandie has been infected by land from west or from the sea. The common C-wave front then progresses significantly to northwest, while the progression to east practically freezes during several months. With small differences, the spatiotemporal calculations closely reproduce the historic evolution indicated by the dates and positions of infected places. Although C-results have a large dependence of [13]-data in the calculations, some discrepancies exist that would deserve discussion. The infection dates of Lyon and places at south-east Auvergne, at April 1348 according to historic data, look like too early, if compared with the surrounding places, reported to be infected later. Numerical results suggest the need for analysis of historical data and model parameters. Paray-le-Monial has calculated infection date at august of 1348, according to a regular spreading of the epidemic. The attributed [13]-date April 1349, quoted as having large uncertainty, lays very much outside the dates of its neighbours. This large difference is shown by a detailed inspection of Fig. 6 and may be related to its position out of Rhône-valley or to the special character of the place, populated mostly by monks at the epoch. The infection of Rouen is mentioned to occur at the beginning of the summer of 1348. If we accept this report by [13], it seems natural Paris and nearby places to be infected two months later. This is compatible with the infective wave to propagate southwards (cf. Fig. 3b), although one expedition in a navigable river at that distance could travel quite faster. In its eastward propagation, the epidemic was frozen in a line that separates cities infected before September 1348 and infected after April 1349, the delay being due to seasons of the year, to political divisions or both. At Alsace and Lorraine, another barrier separates the cities infected before August 1349 and after November 1349. A greater problem arises relative to the quick infection of Rouen from all possible sources, located at 1000 km, like Bordeaux, or even further. Bordeaux is referred as having been infected in August 1348. Sea vessels could normally travel, at that epoch, at an average of 40-50 km/day, even if good conditions could result in greater velocities. The greater the sea travel times were, the greater the probability of arriving nowhere, as consequence of the effects of the disease in the functionality of the crew. From all the possible explanations for the hypothesis of Rouen being infected from Bordeaux, a simple one is the possibility that Bordeaux has been infected via Garonne River at an earlier date than the one mentioned by [13]. The city of Agen, some 120km up by the plane Garonne River, was infected at May 1348, coherently to other dates around. This could allow Bordeaux to be infected before August, making it possible Bordeaux to be the starting point of a quick arrival of the plague to Rouen and then to Paris, these ones possibly at dates not very different from the ones cited at [13]data. C-results with the chosen parameters did not succeed in having one sooner arrival to Bordeaux than [13]-data did. Apparently only locally altering the parameters of equation 3 would it become easy to rise the wave speed. The alteration of D values in the trade route to the city did also not solve the problem, what sometimes happen due to the effect of diffusion perpendicular to route direction. If the hypothesis of one fast travel by the river is correct, the reintroduction in the model of convective transport of large speed would with certainty work. Another possibility involves the introduction of the 1D domain referred in section 2.7 that can preserve the (in some extent artificial) diffusive approximation with smaller transverse fluxes. Although [13] state otherwise, the beginning of June as infection date of Bordeaux would be compatible with the references of [13] about the trip of princess Joan from London, and the possibility of her infection in Bordeaux. If princess Joan left London at mid-July, hardly a boat leaving Bordeaux at the end of June would carry the bad news to the royalty. Moreover, the last 15

vessel from Bordeaux arriving at London before Joan departure could possibly leave Bordeaux some time before the end of June, depending on the dates at which those voyages took place. The royalty may have not be informed of the facts, in the case that Bordeaux infection was not obvious at the time of ship departure. So Bordeaux infection could have started at the beginning of June or even before, hypothesis that have to be tested and confirmed or denied by historical analysis. The reunion of the wave fronts from south and from north could occur, as the historic data and numerical results suggest, after October 1348, near Orléans and Auxerre, according to simulation results and compatible with [13] analysis. Between November 1348 and April 1349, the eastern epidemic wave front does not progress from France to Germany, at a line near the frontier of France with the Holy Roman Empire. Also the progression of the front from Italy to the north freezes during that time. The presence of political and geographic barriers may be responsible by this delay. These barriers are included in the lines of low diffusion coefficient of Fig. 1. The time interval also corresponds to the 1348-49 winter, their effects possibly adding. The existing data about Givry mortality ([13] and references there) provide a very detailed knowledge of daily death rate. The comparison of our numerical results (Fig. 4) with the ones shown in [13] reveal one similar general shape of the mortality rate curves. Nevertheless, the important differences between our results and those data are clear in respect to peak mortality, dates of mortality beginning and the duration of significant mortality rates. The clear [13]-data show a mortality rise at end of July 1348, finishing at November of the same year. Our results indicate the mortality rise at about January 1349 and the end 8 months later, grossly doubling the duration established in [13]-data.

These results suggest that our model may possibly be appropriate to disperse populations but not to high population concentrations. The population density in the region of Givry used in our simulations is about 17 inhabitants per square kilometre, while it is expected that the density inside the real city can be larger than hundreds per km2. Of course if one S-independent contagion coefficient instead the one of equation 3 has been used, the expected results would be highly reversed. In this case, and using local susceptible density of the epoch, one much quicker mortality evolution would occur. The expectation of the present author is that, if compatibility is possible to achieve in one model of this type, major attention must be directed to the mathematical forms of the contagion rates and the values of the parameters. In the most benign hypothesis, the change of parameter values inside equation 3 will lead to realistic results both inside cities and at countryside. These changes may possibly have effects in mortality coefficient, which in turn involve questioning the disease type or types and respective rates of pathogen transfer. This appears to be an important limitation of this model provided with the chosen parameters, that must be further investigated, what is not done here. Nevertheless there is the possibility that the double or triple subdomain option referred in section 2.7 can solve the problem with no need of abandoning this type of models.

Iberia

16

At April 1348, Black Death go beyond the Oriental Pyrenees border and propagates south by the coast, possibly coinciding with the infection of ports from sea. At summer, the north Pyrenees infective wave passes the west Pyrenees edge. By this time, the westwards wave progressing in south Pyrenees approaches the bay of Biscay. Both wave fronts unite after November 1348 and propagate to west. Meanwhile, at the beginning of summer, Galicia receives the epidemic from sea, which then travels to south and to east. Before the end of the year, [13]-data indicate the formation of a broad band by the north of Iberia Peninsula, when the western front joins the wave coming from the Pyrenees. Possibly Palencia like Zamora were infected by the north-western front. In [13]-interpretation apparently these two cities were immersed in the eastern front coming from the Pyrenees. These two interpretations need some more historical work to ascertain the infection date of Palencia and locations eastwards, to which the results from the present model must in principle be subjected. Apparently, a southwards Galicia wave came to Portugal while an infection wave from Lisbon travels north, the two fronts uniting at the end of 1348. Apparently with origin in Lisbon, the epidemics moves south and to Alentejo and Spanish Extremadura. With the parameter fields used, the calculated infection of Coimbra could have its origin in the front from Galicia. The infection wave from Lisbon or from the seaport of Figueira da Foz are also strong possibilities, which only more detailed historical data can decide. The infection of Évora could have occurred around October 1348, from the probable origin at Lisbon, according to the simulations. This date disagrees with the year of 1349 of [13]. Nevertheless the scarcity of data does not guarantee the correctness of the calculated date. The infection of Silves could also have been caused by the same wave front at dates around march 1349, or by an independent wave with origin in Algarve seashore. The scarce historic data do not suggest, neither the numerical simulations do, Silves to be infected by land by a wave coming from Spain. The fitting of C-results to [13]-data are many times very dependent of interpolations or extrapolations in areas with scarce information density, more historical data being required. The wave front of the Oriental Pyrenees progressed to south. Probably, together with the infections of seaports, the Iberian east coast wave arrives at south Spain at the beginning of summer 1348. Portuguese and Spanish southern fronts merged at the summer. The kingdom of Granada could have been infected by sea, according to [13]. The particular simulation realization, with its own parameters and particular conditions, also allows Granada to have been infected by land. Almeria was infected at May 1348. Military confrontations at Huelma in 1348 between Toledo and Granada forces could have led to the transport of the disease agent to Cordoba and then to Toledo. Cadiz, infected at unknown date, cannot also be discarded as the possible infection source. The infection of Toledo is determined by one C-interpretation in a way to which not much credit must apparently be given. Toledo stays at the centre of one definable area of about 200000 km2 with no sound [13]-dates excepting Toledo itself. Relatively small differences in the choice of initial densities and on diffusion parameters in the defined routes and barriers may lead to substantially different results. Ireland

The eastern side of Ireland was soon infected since august 1348, according to historical sources, [13,20]. This early date raises the question of knowing the external source of the transported pathogen to the island. Aquitaine, Galicia, Normandie, Bretagne are possible origins of Ireland 17

infection. Although some historical data point to infections in locations of Great Britain as soon as June 1348, as in the cases of Weymouth and Bristol, in the present study these dates were not considered, as these cities appear to be surrounded by a significant number of very near places infected later. Despite the relatively short distance, it took about one year to infect, possibly by land, the western part of Ireland. Accordingly, a diffusion barrier (Fig. 1) was introduced in the model parameterization to delay the calculated infective wave. If geographic, political or other kind of barrier really existed has not been analysed. Switzerland

[13]-data shows the infection of south Alps generically occurring before December 1348, while places at north being infected after May 1349. The present model incorporates a diffusive barrier with low D values located in the north slope of Alps. Apparently the pathogen enters Switzerland by Geneva and Nyon coming from Lyon. A couple of months later, Bellinzona is infected by the wave passing Varese. The contagion of the centre of the country seems to have occurred with origin in SW and SE, while the contagion of the north can either have origin in the centre of the country or from the west border with France. The complexity of Alps orography suggests the need of markedly more detailed data about initial population density, diffusive routes and barriers. The level of the details used in Switzerland is similar to the level used for the rest of Europe. The consequences are significant doubts about the precision of Swiss simulation results. Specific historical details about time dependence of the transport terms, related to travels and obstacles, could also be needed to arrive to more accurate results. The analysis by [13] leads to more detailed description of the pest evolution in Switzerland. Lucerne and Disentis C-dates are in perfect agreement with [13]-data, while Engelberg is hit by the calculated epidemic wave at a date quite different of the indications of [13]-data, this being supposed to be caused by the limitations of the model details. Austria

Southern Austria was infected from northeast Italy at the autumn of 1348. The extreme northwestern part of the country was infected later, probably by the wave front coming from Switzerland and south-western Austria. The northeast could have been infected from any direction, with the probable exception of north. Waves from southern Austria, from the eastern Germanic border, or even by one wave, not sufficiently reported by historic data, coming from southeast, could have been responsible by north-eastern Austria infection. Neuberg is designated infected at September 1348. At this time the calculated wave is yet far from the city location, Fig. 3b, which only came to be reached at December 1348. Comparing with other nearby places, [13]-data appears to indicate one too early infection, unless could it have been infected from SE, about what there is no information in the accessed data. Vienna received the pathogen at April 1349. C-results indicate May, what seems reasonable when compared to dates of nearby cities and agrees with the presence of the hypothetical diffusion barrier of Fig. 1.

Great-Britain

18

Weymouth port was infected, according to the option of [13]-data, at August 1348. A large number of places around and very close to Weymouth were infected later. It seems likely that the infection of the city occurred at a later time, so September 1348 was included in the calculations. The delayed infection time for Weymouth introduced in calculations, was already used for Venetia. This procedure results in much more regular spreading of the epidemic, and better agreement between numerical results and the historic data. The existence of one delay in the contagion dates reported by historical data is suggested by [10]. Corrections applied to some coastal cities made in the present simulations is applied in just a few cases. From southwest the epidemic spread along the south coast, progressing then to north mainly by the east coast and, a few months later, also by the west coast in its way to north. The infection data for Wales only mention infected places by the coast. Historical data for north Scotland are also limited. Norwich infection in January 1349 is inserted in the calculations as an independent source, presumed arrived from the sea and not as being hit by the interior epidemic wave. London was, according to historic data, possibly infected in September 1348. Historically originated data do not show any infection of London’s neighbour cities before November 1348. The progression of Black Death that these data show seems to come from west. Byrne [10] claims London to be infected at late summer or fall of 1348 via the road from Oxford. In their work, [13] mention the existence of several concerns about time accuracy. The present model do not set the infection of London as an independent sea infection, rather calculating it as a normal land infection. If sea originated infection happens, it is probable to have occurred at times close to land infection of that region. The numerical results indicate the city could have been reached by the epidemic wave at January 1349. A barrier from Newport to the east coast near London seems to exist, only being crossed by the end of 1348-49 winter. A low diffusion line is defined as shown in Fig. 1. More detailed studies may clarify the importance of this particular winter and of barriers related to geography or to political and administration rules.

Scandinavia

According to historical data, Oslo region was infected from April 1349 on, progressing significantly towards NW and SW only after the end of the summer. The reasons associated to these facts or results can be diverse. If [13]-data are accurate about these details, one different time dependence of D on time could perhaps be assumed replacing equation 15. From Norway, the near regions in the SW of Sweden would have been infected. Apparently, the infection of Halmstad could possibly infect the remaining Swedish coast till Uplandia. Places with historically uncertain infection dates like Orebro, Stockholm and Södermanland region could have been infected, according to the numerical results, after June 1350. The calculated epidemic extinguished in part of the peninsula due to the low population density.

19

Belgium

[13]-data strongly suggest the appearance of the Black Death in Belgium from nearby French cities. These cities were reached by the coastal western wave to which the wave with focus in Calais united (Fig. 3c). Between France and Belgium, one 100 km barrier was imposed to deal with the time delay of infection between Calais at m=15, Dec48, and nearby Belgium cities, infected after m=20. The progression in Belgium occurs regularly to east and south. This wave may have infected Holland, but this hypothesis has not enough support from [13]-data.

Germany

One barrier at Alps northern slope, mostly in the interior of Switzerland and Austria, divide southern and northern regions. The north was infected typically 5 months after the infection of the south. After this 5 months delay, the wave front from west Austria infects south-eastern Germany at the beginning of 1349 summer. Moving westwards, by August 1349 this wave front meets at BadenWürttemberg the infective wave coming from Switzerland or from Alsace, France. Apparently, the wave from west advanced at summer beginning, while the southern wave passed the Alps barrier and, in less than two months, they arrived to the barrier inside Germany (see Fig. 1), more than 150km far from the Alps barrier. In this process, the regions of Baden-Württemberg, south Bavaria and the Czech Republic, except northwest, were quickly infected according to [13]data. Typically, all referenced cities in the band between the two barriers were almost simultaneously infected at dates very close to July 1349. This includes Rhineland-Palatinate, Baden-Württemberg and the southern part of Bavaria. This historically reported simultaneity seems quite improbable to take place in a so large geographical area, [13] referring inaccuracies in many data concerning Germany. The simulations performed indicate the interior barrier to be impervious, the epidemic being overtaken by its borders, the epidemic then moving southwards after February 1350, finally reaching the inner angle near Wurzburg and Nuremberg after June of the same year. This go around movement is not denied by [13]-data although, due to the low spatial information density and the low reliance in it, this description do not constitute by now one established hypothesis. In the even vaster area at north of the previous barrier, the reported dates of infection typically lay between February 1350 and August of the same year, showing no clear spatial sequence that could undoubtedly indicate the evolution of the infective waves. There is the possibility that this vast area could have received the plague from different directions, namely from Rhineland-Palatinate at SW, from the north harbours and from the south of the barrier inside Germany, in this hypothesis contradicting C-results. Within the context of insufficient data, Schleswig, near the coast, was imposed as an independent maritime infective source. Czech

The very scarce available data allow the possibilities of Czech infection from south or southwest, not being possible to exclude the hypothesis of contagion from southeast. Netherlands

20

The few available data in principle allows for the possibility of infection by a northern harbour, from the eastern border with Germany, from southern border with Belgium, or even by combinations of these possible sources. In the numerical simulations, the arbitrary option was taken of not including the port of Foswert as one local independent epidemic source, so linking its infection to the wave coming from the western coast. Zwolle and Foswert are infected at [13]-dates differing 11 months although separated by only 70 km. It is probable that inaccuracy in historical data may be responsible by this time gap. Nevertheless, if these or similar dates are confirmed, it will be necessary to find the possible ways by which this time separation is so different from the usual ones in this epidemic. In the case of Swiss cities of Engelberg and Dissentis, one similar delay happens, which can be explained by barriers to infection waves. Inland non-navigable water could be the barrier, if it was there at that time. In Holland case, if dates from [13]-data are correct, the german origin of the infection of South Holland is clearly favoured. In such case, the inclusion of a barrier to the transport of the infective agent should have to be imposed between the two mentioned cities. The delay of 5 months between infection dates of Zwolle and Deventer, separated by only 40 km, would need the further imposition of a diffusive barrier. This level of detail in one region with so few historical data was not included in the model. The fight between regularity and specificity is present in many situations and is not quantified, being here solved in the benefit of regularity. Denmark

From the few historical data available, Zealand Island and Jutland may have been infected in the same month from sea. Jutland could also have been infected by land from south. The results of simulations show Alborg to be affected by the land epidemic wave near the same month indicated in the historical data. This does not prove that the city was not infected by sea. Poland

In the historically based data of [13], Poland was the last European country infected by the epidemic to be studied. Only a small region of north Poland has some historical information about infection dates of cities. Gdansk is used in the simulations as an independent infective source. The supposedly most populated Vistula basin was infected at a relatively high calculated wave velocity, the imposed high diffusion line (Fig. 1) also contributing to this effect. Eastern Europe

The data from [13] utilized in the present study do not include a number of countries in the Eastern Europe, possibly due to deficit of historic sources, though the plague has been present there. With the possible exceptions of Finland and Iceland, [7], there exist data that indicate the epidemic spreading in almost all European regions. The majority of the published historic data does not include enough information for significant parts of Europe, as [27] refer.

In Fig. 3, the maps from [14] show blank areas, that correspond either to uninfected regions or regions about which there is not enough information. The mortality coefficient in Table 1, leads to longer reach of the infected individuals than in the hypothesis that Black Death is a bubonic plague with a mortality coefficient about =0.09. In this 21

latter case, more uninfected areas arise from simulations, but not enough knowledge exists about real blank regions to take conclusions about the epidemic type, based alone on this aspect. In Fig. 5 the overall picture of pest progression is shown for C-results and for [13]-data.

Fig. 5a, simulation results, display an increased regularity of the pest evolution, compared to Fig. 5b that shows details that are not discriminated in C-results. 3.2 Error assessment One essential guide to the parameter optimization has been the measure of the fitting quality between simulation results and [13]-data. Comparison of the infection [13]-dates of cities with the homologous C-results allowed a laborious choice of parameter values to diminish errors. Automatic calculations by Monte Carlo methods to optimize parameters for better fitting revealed problematic. An enormous number of data would have to be optimized in order to reduce the error. Besides the parameters of table 1, also the criteria for epidemic existence, for epidemic declaration and the initial population density would have to be optimized. The choice of the boundaries of parameters values, the limits to local population density in order to not in large scale contradict the available values, are important parts of the problems that have been faced. Of course it is possible to model the parameters limitations to minimize errors and retain some criteria. But it can be an excessive work for the reward. The utilized [13]-data are taken as the reference values to comparisons and error calculation. Meanwhile, as referred in [13], the accuracy of many of the direct historical data is very low. Even using data from miscellaneous sources, their quantity and quality may be scarce and limit the confidence in some of the final results, even after the use of methods that in some way filter, blend and improve data to obtain the final product. The comparison of C-results with [13]-data has to be considered as one verification of the pair mutual agreement, although in practice it will be mostly used as the verification of the quality of Cresults in relation to [13]-data. The main parameter used to evaluate error will be the mean absolute error  of the dates of infection at cities locations. The standard deviation is not preferred as it gives relatively large weight to large deviations not resulting from a very large number of independent factors. Cities with small precision in the historic infection date were not included in the error calculation as in the cases of Évora and Badajoz, for which only the year is specified. Lund, declared infected at the summer of 1350, was included in the calculations assuming infection at July of the same year. At November 1351 (m=50), 363 cities having been chosen from [13], the mean absolute error between dates of cities infection is =1.54 months. The regularity that models like the present one normally project in the results, what also happens in the present study, and the relatively small differences between C-results and [13]-data, suggest the idea that also [13]-data can be considered regular in the sense of showing close to steady progression of the epidemic over wide areas. The similarity between C-results and [13]-data is however no proof of the quality of both data, as a substantial amount of initial C-data is extracted from [13]-data and parameter optimizations are performed just lo reduce the error. It is sure, in the absence of more information, that, in many cases, the option to choose the irregularities of [13]-data or the regularities of C-results is an uncertain job. Nevertheless, more credit is in principle attributed to [13]-results than to C-results, despite surely existing some cases in which the opposite must be correct. 22

In Fig. 6 the advances and the delays of simulated dates of infection in relation to [13]-dates of cities is mapped. For most of the area, Fig. 6 shows the difference between dates of arrival to be within -2 and 2 months, where negative numbers indicate that C-spreading arrive sooner than [13]-dates of city infection. C-infection is determined under the criteria of section 2.6 for epidemic perception. Large positive or negative shifts in (tC-tH) happen at some cities. The map of Fig.6 can be viewed as a measure of the error of C-results as compared to [13]-data. In reverse, it can be used as an estimator of the irregularity of [13]-data in respect to the dates of pest beginning. Both views may serve to seek reasons for the irregularities in (tC-tH). 3.3 Number of cities infected at each month Fig. 7 shows the total number of places infected till a certain month from [13]-data and the number of places that have already been reached by the calculated epidemic waves at the same month. The overall agreement is clear, even if for a large number of months there are more places reached by the calculated wave, while for other months there are more cities infected as indicated by the historical data. Fig. 8 shows the number of infected cities in each month, according to historic data and from calculated results. The different local progress velocities between historical and simulation results is responsible for the advances and delays of the instantaneous maxima and minima. The overall differences between the behaviours of historic data and calculated results are nevertheless relatively small as shown by Fig. 7. The overall trends of the infection rates of the number of cities appear to indicate more infected cities from spring to autumn. Some care must nevertheless be taken. The observed maximum and minimum rates at October 1348 and February 1348 seem to deserve some research. The absolute maximum number of 28 infected cities in one month is observed at m=13, Oct48, and is followed by one steep diminishing till 2 infected cities at m=17, Feb48. So, from the very high value at m=13, almost linearly goes to one of the smaller levels of all the time of Black Death, rising again suddenly to 23, in just one month, one of the highest levels at m=18, Mar48. At m=16 and 17 the epidemic waves are mostly located at the great continental barrier and at the GB barrier, at times when progression is very slow. So very few cities were infected at those months. In Iberian peninsula, the wave fronts are located at one vast region where the density of cities with identified infection date is very low. Also in this case the number of new infected cities with known [13]-dates is very low. The huge amount of infected cities at month 18 does not lay at the north and east sides of the great continental barrier. This large amount is located at GB in the regions of East England, West Midlands and the south part of East Midlands, where more than half of all European cities with identified infection date at March 1349 are located. In one month the distance the wave traverses is more than 100km, in one area with a large density of cities with known infection dates. And the absolute maximum city infection rate at month 13 is also peculiarly distributed between wave fronts at Spain and south England, a few in the rest of Europe. 3.4 Sensibility of error to parameters

23

The model development had in mind the main dynamic features of the epidemic under study. Initial conditions and the values of the parameters were subject to prior estimates and then adjusted to obtain better fit to the cited [13]-data. The reliability of the results depends on the model, initial conditions, coefficients and some utilized criteria related to infected density. All these components may be improved, preferably in an interdisciplinary approach. Although in principle possible to all items, the calculations of the effects of variations in the final errors will be done relatively to just some parameters. The most important parameters used, Table 1, determine the final results, which are sensible to the set of chosen parameter values. Important errors are present in the original historical sources that are not confined to small intervals, namely about the dates of pest local start. Although filtered by the methods applied in [13], the final [13]-data inevitably carry uncertainties and errors. These data will be used as reference for comparisons owing to the condition of apparently being the more reliable set of organized quantitative data on this subject. Figure 9 shows the calculated effects of parameter values in the average absolute error between start dates of infection of cities obtained by the present study and in [13]-data. The common position of all minima means that approximately the best values of these parameters were chosen, other initial conditions being fixed.

In some cases, small modifications of parameter values can lead to some changes in the evolution of calculated waves, leading a small number of cities near wave fronts to be included in or excluded from the calculated infected area. This may originate some noticeable error variations even for small modifications of parameter values, which is the case for kc0 variation in Fig.9. 3.5 Distribution of the arrival time errors For each individual city, C-results and [13]-data frequently differ. The calculated date of infection at one city location happens m=mC-mH months after the prescribed by [13]-data (m<0 if C comes earlier). Which cities have m advances or delays is implicit in Fig. 6, if it has enough resolution and if we know the location of each city. The number of cities that have mC-mH equal to a certain number m is represented in Fig. 10 for all differences m found from C-results and [13]-data. For comparison, the corresponding to the Gaussian distribution is shown in Fig.10, normalized to sum the total number of cities. In the cases of infection of harbours from abroad, in which mC has been imposed to be identical to mH, this type of cities was removed from the data used in the figure. The values of average m and standard deviation were calculated on the basis of the data relating C-results and data from [13]. The Gauss distribution curve in m space obtained using these parameters has been normalized so as the total areas (meaning total number of infected cities in question) to be equal. The parameters that define the Gaussian curve are the average value of m=-7.95x10-2 and standard deviation =2.26.

The lower peak value for the Gaussian function and the shift between peak locations suffer the influence from the tail at low m where the calculated dates for a small number of cities were much sooner infected in relation to the dates in [13].

24

While the absolute mean error of mC referenced to mH takes a value close to 1.5 months, the measure of the error by means of the standard deviation throws the error to a value larger than 2 months, due to the influence of 2 or 3 cities isolated in terms of m. The importance of large m deviations, like in the case of Zwolle with a difference to the average of 12 months and two other cities with deviations of 8 months between mC and mH are clearly located outside the bulk of results. The option of ignoring those cities was not taken, although they are clearly outside the bulk of date in m scale. These deviations from the normal distribution led to the option of considering the absolute mean error as the measure of error to be adopted here. 3.6 Mortality maps Fig.11 shows the maps of monthly mortality rate obtained from our simulations. The comparison with the ones obtained by [13] for cities shows one considerable delay in the evolution of our calculated results and a greater mortality rate. The delay in the C-results can be related to the smaller density in the countryside. In the present model all places are equivalent to countryside. The delay of infection of countryside relative to the cities located at trade routes is expected to be a real phenomenon that is not represented in the present model, although referred in section 2.7. Also specific criteria for infection beginning at cities and countryside may be involved. Low mortality coefficient of the infected, resulting from large time from infection till death, allows a longer period for contagion to occur. In such cases, high initial susceptible densities lead to large mortality for long times. If this large total death fractions really exist in Black Death, the hypothesis of a relatively low mortality coefficient seems to be more probable than one high coefficient. In the opposite case, if large mortality in the more populated regions does not happen, a quicker death, larger mortality coefficient, may be the cause. These several questions must be addressed in trying to solve the discrepancies of Fig. 11. 3.7 Correlations between mortality and infection Some graphical representations linking infected density and death rates can be useful to predict possible evolutions of the epidemic and expected time delays. Fig. 12 shows the lines of initial perception of pest at certain months (according to the particular criterion of section 2.6) and the lines of maximum mortality at the same or other dates. As expected the lines of pest perception come much earlier than the ones of maximum mortality. At some places, like in Balkans, the perception of pest beginning at month m=10 in some places coincide with the places of maximum mortality at the later month m=20. In other words, at certain places the maximum mortality happened 10 months after infection perception, according to this mean-field model with the inserted parameters. Looking simultaneously at some different places at month m=25, maximum mortality already happens at south Adriatic, while central Greece, located 300km SE and belonging to the same calculated propagation wave, is beginning to have the perception of the regional tragedy that is going to develop. These results are obtained under the conditions imposed to the simulations, that presume the infection came from north in this region.

Other quantitative links between infection and death rate can in principle be built, based on equation 12. The rate of contagion is related to mortality, although with the presence of other influences, like the existence of infected fluxes and the changes of susceptible density. The comparison of month mortality with kc.SI may reveal to be more clarifying than directly with values of I. The behaviour of 25

temporal and spatiotemporal correlations between terms like these can help to describe some particularities of the epidemic and serve to make quick anticipations. The relation between month mortality fraction and infected density at different times is shown in Fig. 13 at one place chosen to coincide with the position of Givry, France. In the present study the behaviour of city locations is not different from the ones at countryside, as no spatial heterogeneity in initial susceptible density is being imposed. Monthly mortality fraction and density of infected begin with close to zero values, at that location, the curves showing the evolution of both values in anti-clock-wise direction while time grows. The relation between month mortality and infected density 24 days (supposed survival time of the infected person) before the calculated month mortality is shown in Fig. 13a. The relation between month mortality and infected density at the end of the same month is shown in Fig. 13b. For comparison, see related results in [13]. Fig. 13a shows the values of monthly mortality fraction to be close to proportional to the values that density I assumes 24 days before. At Fig. 13b, with mortality and I considered at the same times, the number of infected, or of its density, at early dates do not immediately result in deaths, that are going to grow in the next months. After the maximum value of I, mortality still grows until both come to decrease till local extinction of pest.

4. Conclusions One deterministic continuous model of the epidemic dynamics was built on the basis of the balances of the epidemic species, with terms of creation, destruction and transport, calculated by means of previously existing and by new submodels to determine the rates and fluxes of the balance equations. The model of the epidemic dynamics is similar to the one of [23]. One term for the transport of susceptible was not included, assuming the phenomenon it describes has not strong manifestation in the epidemic dynamics. The transport of infected was found important at trade routes and was included in a diffusion form in the appropriate balance equation. A new term relating contagion and pendular short range dislocations of individuals was included. The rates of infection resulting from the temporary displacements of susceptible and infected between locations with different densities of susceptible and infected are calculated. These terms of interaction of the pendular motion with contagion have diffusion-like forms. The contact rates between humans show some saturation effects that influence the contagion rate. To represent the phenomenon of contagion saturation, modifications to the classical expression of the contagion rate has been introduced. The modifications preserve the proportional dependence of the contagion rate with the infected density, modifying the dependence on susceptible density so as to achieve the saturation effect. One term of mortality of the susceptible, active after the epidemic wave front arrival, is included to represent the effects of economic and social crisis on death rates. Most of the parameters included in the equations have constant and uniform values. Some of them are modulated, like the diffusion coefficients, which have greater values at summer than at winter. The diffusion coefficients are increased at the lines roughly coinciding with the trade routes, and are diminished at lines where geographical or political barriers exist.

26

The type of pathogenic agent or agents responsible for disease are important in the effects in mortality coefficient. The option was made, after [14] and references there, to not consider Black Death as one bubonic plague, with effects on the value of mortality coefficient. The initial susceptible population density was taken from demographic references and modified in a way to enhance the agreement between the simulated results and the historical data. At specified locations and times, infection sources were imposed, modelling the effects of the arrival of the pathogen by long range sea travels. Threshold values for the values of infected density were imposed, below which the epidemic becomes locally extinct. One place is considered infected when the infected density exceeds one value associated to the local perception of epidemic existence. This declaration is used to compare the dates of city infections calculated by the present model with the dates determined at [13]. The differential balance equations have been approximated by algebraic equations by means of the finite volume method. With the initial and special conditions, the algebraic equations, supplied with constant or variable parameters, were used to simulate the evolution of the epidemic system as described by the model. The results obtained on these basis were compared with the results of multidisciplinary origin of [13]. The comparison of simulation results with organized external data revealed the capability of the model to reproduce some of the main features of the spatial-temporal epidemic spreading. However, in areas with limited historical information, the simulation results cannot be adequately compared with historical data, what happens especially in east and southeast Europe. In some regions, like Switzerland, the initial conditions and the parameters would have to be much more detailed to arrive to reasonable calculated results. The spatiotemporal evolution of the calculated epidemic spreading closely reproduces the epidemic evolution reconstructed from [13] data. The calculated dates of infection of cities were compared with the historic homologous, the mean absolute error having the value 1.54 months. This low value shows that the present mathematical model has the capability of making simulations whose results are in close agreement with some of the available data. The model of the contagion–pendular motion interaction appears to adequately reproduce the phenomenon, partly avoiding the, in some cases insufficiently sustained, classical gradient diffusion term. The calculated global mortality of 61% is outside the majority of historical estimations, typically between 25% and 50%, although some estimations exist with higher values. Simulations using mortality coefficient corresponding to bubonic plague led to results close to 33% of mortality. No conclusion or suggestion is extracted from these results about the nature of the involved pathogens. The adopted mortality coefficient corresponds to diseases with an average survival of about 24 days. One higher , would in principle lead to lower global mortality and give rise to more uninfected areas. The results from Christakos et al [13] present areas with no epidemic or with uncertain presence of epidemic. The modifications of the coefficients could then lead to results with regions free of pest that are nevertheless not indisputably free of pest. This presence of uninfected regions in Western Europe happens in the case of simulations under the assumption of Black Death to be one bubonic plague. The large uncertainties remaining about the model results do not suggest any option about the diseases acting in the epidemic. 27

Local time profiles of calculated mortality, surely also of infected density, are different from the existing profiles for cities. This feature is supposed to be mostly related to the dynamical differences between low quasi-uniform densities defined in the model and the high density inside cities. Also definitions of parameter values could influence. Alterations suggested in section 2.7 to the model can in principle modify the results to one yet unknown extent. These alterations would allow the simulations to cope simultaneously with countryside, trade routes and cities, inserted in one unique model with three interrelated spatial subdomains. Some groups of values attributed to the parameters and initial population density, somewhat different from the chosen ones, have been tested, leading to similar results of the velocity patterns of the epidemic waves and total mortality. The choice of well-founded quantification for all the parameters, as well as the initial susceptible density field, must be undertaken by means of interdisciplinary work. The ensemble formed by the model, the field of initial population density, the parameter values and the initial infection places at specific times, was found to have the capability to simulate the evolution and effects of the Black Death, with some results in good agreement with the historical data. A wider ensemble also including historic data can be analysed. This type of work has already been developed by [13] using different methods. The inclusion of models like the present one can also be done, possibly leading to combined improvement of some of the ensemble elements. The possibility to further execute detailed evaluations of the quality of each of the elements of this ensemble, their combined improvement and hypothetical stabilization is far beyond the capabilities of the present study. Wheeler’s ‘R’ may suggest that reality is built by the observer. It is surely better to maintain the concept of the existence of an independent reality. In this sense, Wheeler’s ‘R’ must be understood as a model the reality we build at every day. If detailed enough we would possibly see the posts to be composed by posts, papier-mâché, much work and something called reality.

References [1] W. Abel, Désertions rurales: bilan de la recherche allemande. In: Villages Desertes et Histoire Économique, XIe – XVIIIe Siècle, SEVPEN, Paris, 1965. [2] T. Andersen, P. Jensen, C. Skovsgaard, The heavy plough and the agricultural revolution in medieval Europe, discussion paper 6-2013, Univ. of Southern Denmark, (2013). [3] V. Artemov, S.B. Beale, G. de Vahl Davis, M.P. Escudier, N. Fueyo, B.E. Launder, E. Leonardi, M.R. Malin, W.J. Minkowycz, S.V. Patankar, A. Pollard, W. Rodi, A. Runchal, S.P. Vanka, A tribute to D.B. Spalding and his contributions in science and engineering, Int. J. Heat Mass Transfer, 52, 2009, 3884-3905. [4] A.W. Bacot, Further notes on the mechanism of the transmission of plague by fleas, J. Hygiene 14 (1915) 774-776. [5] P. Bairoch, J. Batou, P. Chèvre, La population des villes européennes de 800 à 1850, Droz, Genève, 1988. [6] D. Balcan, B. Gonçalves, H. Hu, J.J. Ramasco, V. Colizza, A. Vespignani, Modeling the spatial spread of infectious diseases: The GLobal Epidemic and Mobility computational model, J. Computational Science, 1 (2010) 132-145. [7] L. Börner, B. Severgnini, Epidemic trade, Freie Universität Berlin, discussion paper 12, (2011). [8] O.J. Benedictow, The Black Death 1346-1353: The Complete History, Boydell Press, Woodbridge, 2004.

28

[9] M.G. Buhnerkempe, R.J. Eisen, B. Goodell, K.L. Gage, M.F. Antolin, C.T. Webb, Transmission shifts underlie variability in population responses to Yersinia pestis infection, PLoS One 6 (2011) e22498. [10] J.P. Byrne, Encyclopedia of the Black Death, CA: ABC-CLIO, Santa Barbara, 2012. [11] B. M. S. Campbell and L. Barry, The population geography of Great Britain c.1290: a provisional reconstruction, in: C. Briggs, P. Kitson, S. Thompson (Eds.), Population, welfare and economic change in Britain, c.1290-1834, Boydel, Woodbridge, 2014, pp. 43-78. [12] S.K. Cohn, Epidemiology of the Black Death and successive waves of plague, Medical History, 52 (2008) 74-100. [13] G. Christakos, R.A. Olea, M.L. Serre, H.-L. Yu, L.-L. Wang, Interdisciplinary Public Health Reasoning and Epidemic Modelling: The Case of Black Death, Springer, New York, 2005. [14] G. Christakos, R. Olea, H. Yu, Recent results on the spatiotemporal modelling and comparative analysis of Black Death and bubonic plague epidemics, Public health 121 (2007) 700-720. [15] J. Edwards, The transport system of medieval England and Wales – A Geographical Synthesis, Unpublished PhD thesis, Univ. of Salford, 1987. [16] A. Hufthammer, L. Walloe, Rats cannot have been intermediate hosts for Yersinia pestis during medieval plague epidemics in northern Europe, J. Archaeological Science, 40 (2013) 1752-1759. [17] J. Kaplan, K. Krumhardt, Zimmermann, The prehistoric and preindustrial deforestation of Europe, Quaternary Science Reviews 28 (2009) 3016-3034. [18] P. Y. Karimi, Conservation naturelle de la peste dans le sol, Bulletin de la societé de pathologie exotique 56 (1963) 1183-1186. [19] M.J. Keeling, C.A. Gilligan, Bubonic plague: a metapopulation model of a zoonosis, Proc. R. Soc. Lond. B 267 (2000) 2219-2230. [20] M. Kelly, A History of the Black Death in Ireland, Tempus, Stroud, Gloucestershire, 2001. [21] H.H. Mollaret, Conservation experimentale de la peste dans le sol, Bulletin de la societé de pathologie exotique 56 (1963) 1168-1182. [22] J. D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1989. [23] J. V. Noble, Geographic and temporal development of plagues, Nature 250 (1974) 726-729. [24] N. Pounds, C. Roome, Population density in fifteenth century France and The Low Countries, Annals Assoc Am. Geographers 61 (1971) 116-130. [25] S. Scott, C. Duncan, Return of the Black Death—The World’s Greatest Serial Killer, John Wiley & Sons, Chichester, 2004. [26] P.L. Simond, La propagation de la peste, Annales de l’Institut Pasteur, 12 (1898) 625-687. [27] C. Tsiamis, E. Poulakou-Rebelakou, A. Tsakris, E. Petridou, Epidemic waves of the Black Death in the Byzantine Empire (1347-1353 AD), Le Infezioni in Medicina, 3 (2011) 193-201. [28] H. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics, Longman, London, 1995. [29] A. Yersin, La peste bubonic à Hong Kong, Ann Inst Pasteur, 8 (1894) 662-667. [30] J. Warnatz, U. Maas, R. W. Dibble, Combustion, Springer, Berlin, 1996. [31] M. Wheelis, Biological Warfare at the 1346 siege of Caffa, Emerging Infectious Diseases 8 (2002) 971975.

29

Biography

Orlando Silva is a researcher of the Physics department of the University of Évora (Portugal). Areas of interest include transport phenomena, either inert or reactive, in continuous media. Systems concerning turbulent reactive flows, river dynamics, earth mantle convection, epidemics, have been studied. For most of them, models have been developed, the simulations of which being done with the use of computational fluid dynamics methods.

30

Fig. 1. a) Lines of high D thought to roughly reproduce the main trade routes in Europe; b) Main

geographical and political barriers with low D values.

Fig. 2. Adopted European population density at October 1347.

31

Fig. 3. a) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at February 1348; European infected area, from [14] (continues).

Fig. 3. b) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at September 1348; European infected area, from [14] (continues).

32

Fig. 3. c) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at February 1349; European infected area, from [14] (continues).

Fig. 3. d) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at September 1349; European infected area, from [14] (continues).

33

Fig. 3. e) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at February 1350; European infected area, from [14] (continues).

Fig. 3. f) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at September 1350; European infected area, from [14] (continues).

34

Fig. 3. g) Calculated infected areas and cities with known or unknown infection dates reached by the infective waves at November 1351 (concludes).

Fig. 4. Monthly mortality fraction at Givry location obtained from simulation results. 0.30 0.25 0.20 0.15 0.10 0.05 0.00 10

15

20

25

30

35

Fig. 5. Maps showing the locations and times (m=1 at Oct1348) at which infection begins: a) results of simulation using the present model; b) based on data from [13].

Fig. 6. Map of the difference between calculated month of cities infection and the dates from [13].

36

Fig. 7. Total number of infected cities in each month, [13]-data and C-results (month 1 = oct47).

400 350 300 hist

250

calc

200 150 100 50 0 1

7

13

19

25

31

37

Fig. 8. Number of cities infected monthly, [13]-data and numerical results (month 1 = oct47). 30 25 hist calc

20 15 10 5 0 1

7

13

19

25

31

37

Fig. 9. Sensibility of the average absolute error of calculations relative to the normalized values of

some used parameters. 1.58 D sigma kc0 alfa miu

1.575 1.57 1.565 1.56 1.555 1.55 1.545 1.54 1.535 0.92

0.94

0.96

0.98

1

1.02

1.04

1.06

1.08

37

Fig. 10. Number of cities for which the simulated infection arrived at dates separated by m=mC-mH

relatively to the historically-based infection dates and the corresponding Gauss distribution. 80

ncities 70

mC-mH

60

gauss

50 40 30 20 10 0 -15

-10

-5

0

5

10

m

15

38

Fig. 11. Maps of month mortality fraction from model simulations and from [14] at different months.

39

Fig. 12. Map of the locations where maximum mortality (black) occur at indicated months; lines

where the perception of infection (red) occur at indicated months.

Fig. 13. a) Month mortality vs. the infection density 24 days before the end of the month; b) month

mortality vs. the infection density at the end of the same month. 0.30

monthly mortality fraction

monthly mortality fraction

0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.0

1.0

2.0

3.0

density I at specified times

0.25 0.20 0.15 0.10 0.05 0.00 0.0

0.5

1.0

1.5

2.0

2.5

density I at specified times

40

Fig. 14. “What we call reality consists of a few iron posts of observation between which we fill in by

an elaborate papier-mâché construction of imagination and theory.”, John Archibald Wheeler.

41

Table 1 Coefficient values used in simulations.

kc0(km2.day-1)

(km2)

(day-1)

cr(day-1)

tcr(day)

D0(km2.day-1)

(km2)

9.28x10-4

0.0075

0.041

0.00010

700

2.9

170

42