数理解析研究所講究録
第 828 巻 1993 年 138-149

138

MathematicaI Considerations on Spreading Phenomenonof
Early Farming in Europe
Hiromi SENO
Deparrment ofMathematics, Faculty ofScience, Hiroshima University
Kagamiyama 1-3-1, Higashi-hiroshima, 724 Hiroshima, JAPAN
先史ヨーロッパにおける農耕文化圏拡大に関する数理モデル考察
広島大学理学部

瀬野裕美

A mathematical model of the dispersal of colonies produced by the stochastic migration processdepending on the total population of the group is considered. The expected velocity of the spatial
expanding of the settlement range of colonies is analyzed, utilizing the fractal concept applied tothe pattern of spatial distribution of colonies. The model is used to consider the spreadmg
phenomenon of early farIninng in Europe, with the data of neolithic sites with C-14 dates.

INTRODUCTION
The expanding of the distribution area of some animals has been theoreticallystudied by mathematical models. As for patterns of spatial distribution and the
expanding velocity, some diffusion models have been applied to understand suchphenomena (Ammerman and Cavalli-Sforza, 1984; Martin, 1973; Mosimann andMartin, 1975; Okubo, 1980; Skellam, 1951). Those phenomena considered bydiffusion models should have such a characteristic that the spatial distribution canbe regarded as continuous in space. However, for such phenomena that the spatial
distribution essentially consists of a number of spatially disconnected islands, thatis, colonies, the analysis by the diffusion model has to require some additional
assumptions, and should be regarded as an approximate approach.
In this paper, for the expanding of settlement area consisting of a numberof colonies is presented, a mathematical model of stochastic migration processes is
proposed (Bartlett, 1978). In order to give the relation between the number ofcolonies and the settlement area occupied by them, the fractal concept is
introduced. Analyzing the model, we derive the expected velocity of the
expanding of settlement area. The model is applied to the data of neolithic sites
with C-14 dates, which was used by Ammerman and Cavalli-Sforza (1984) inorder to discuss the spreading phenomenon of early farming in Europe. The
expanding velocity of the settlement area of farming colonies is estimated by ourmodel.

139
COLONY PRODUCTION
Free Migration Process: A new colony is produced by a random migration
process in the existing group of colonies, with a constant migration probability
independent of any other parameter (Bartlett, 1978). That is, the probability of the
production of a new colony is constant, independent of any other parameters.
Now, it is assumed that the colony does not become extinct once it is produced.
Under these assumptions, the following model can be defined:$\frac{d}{dt}P(k, t)=-\lambda P(k, t)+XP(k-1, t)$$P(k, 0)=6_{k)}$

(1)

,

where

. the probability of production of new colony during $(t, t+dt)$$P(k, t)$ : the probability of colony productions during time-period $(0, t)$ .
$\lambda dr$$k$is the kronecker’s delta so that the initial condition means that there is nocolony production at $t=0$ . This colony production system results in the Poisson
probability distribution $P(k, t)$ :$6_{k0}$(2)$P(k, t)= ff^{\lambda t}\frac{(\lambda t)^{k}}{k!}$The expected number of colonies produced during $(0, t)$ is$k h=\sum_{k=0}kP(k, t)=\lambda t$(3),and the expected time of the k-th colony production is$\#k=\int_{0^{\infty}}\tau P(k-1, \tau)\lambda d\tau=\frac{k}{\lambda}$

.

(4)

140

.el.byScheemiaz 線 crfietni 躍註灘欝露盤 ent$rangFig$

explanation, see text.

Size-dependent Migration Process: The migration probability is assumed todepend on the total population of the group (Fig. 1). This means that the colony
production is enhanced more and more as the total population becomes larger.
Under this assumption, we consider the following model:(5)$\frac{d}{dt}P(k, t)=-\}\iota N(t)P(k, t)+\mu N(t)P(k-1, t)$り,$P(k, 0)=$where

: the total population size of the group of colonies at time$\mu N(t)dt$ ; the probability of production of new colony during time-period $(t, t+dt)$ .

$N(t)$$t$This colony production system in time results in the Poisson probability
distribution $P(k, 7)$ in time :
$t$

$T$

(6)

$P(k, T)=e-\downarrow 1T_{\frac{(\mu\tau Y}{k!}}$where the time$T$

is now transformed from time as follows:
$t$

(7)

$T=T(t)= \int_{0^{t}}N(\tau X1\tau.$

Since a colony does not become extinct after its production, we find thatas
. The above result in time coincides with that for the previous case, that
$Tarrow\infty$$tarrow\infty$

$T$

141
is, for the case of free migration process. Therefore, the expected number ofcolonies produced during $(0, t)$ is$\Psi\lambda=\sum_{k=0}^{\infty}kP(k, T)=[\iota T=\mu\int_{0}^{t}N(\tau)d\tau$,(8)

and the expected time of the k-th colony production is$t^{\tau}k=\int_{0}^{\infty}\tau P(k-1, \tau)\}\downarrow d\tau=\frac{k}{\downarrow 1}$.
Then, the expected time in

$t$

(9)

can be obtained through the following relation:$t k=T^{-1}(\{\phi=T^{-1}(\frac{k}{\mu})$

(10)

where

$T^{-1}$denotes the inverse function of $T=T(t)$ .

EXPANSION OF SETTLEMENT AREA
Next, we consider the settlement area of the group of colonies. The settlement

area at time corresponds to the area that has been occupied by those existing$t$colonies at the time. We characterize the settlement range by the minimal
diameter, say , which can include all existing colonies.
In the case when the settlement area expands in every direction with thesame probability, the shape of the settlement area can be approximated by the disc,and therefore, when the spheric nature of the earth can be negligible and be
approximated well by the plane, the range approximately has the following
relation with the total number of colonies $M;M\propto R^{2}$ . However, since the
$R$

$R$

expanding of the settlement area is constrained by geography, climate, cultural
factors, etc., the shape is in general possibly inhomogeneous in direction. It is
likely that the shape has a fractal nature (for the concept of’ffactal’, see, forinstance, Mandelbrot, 1982). To deal with such cases, we assume the generalized
relation between the settlement range and the total number of colonies as follows:142

Fig. 2. Illustrative explanation of the relation of the fractal dimension to the spatial pattern ofcolony distribution. Each black disc shows each colony. (a) $d\sim 1;(b)1<d<2;(c)d\sim 2$ .
$d$$M\propto R^{d}$$(1 \leq d\leq 2)$

(11)

,

where the power charactenizes the spatial pattern of the settlement area occupied
by colonies (Fig. 2). It is called cluster dimension or mass dimension, which isone offractal dimensions. When $d\sim 2$ , the spatial distribution of colonies can beapproximated by a disc. When $d\sim 1$ , the distribution is one dimensional, that is,$d$the colonies are arrayed on a curve. For example, the latter case may correspondto the case when the colonies are located along a river.
Through the relation (11), we can consider the velocity of the expanding of
the settlement range. That is, the velocity is given by$V$$V= \frac{dR}{dt}\propto\frac{d}{dt}(M^{1/d})=\frac{1}{d}\cdot M^{(1-yd}\cdot\frac{dM}{dt}$

.

(12)

Since the expected tota1 number of colonies at time is given by $kh$ , the expected
range of the settlement area is proportional to. Therefore, we consider the
expected velocity at time as follows:$t$$k\Psi^{d}$$t$$V_{t}$$\ovalbox{\tt\small REJECT}\propto\frac{1}{d}\cdot\phi h^{(1-dy_{d}\ovalbox{\tt\small REJECT}}d_{dt}$

.

(13)

For the case of the free migration process, the expected velocity is$\overline{V_{t}}\propto\frac{\lambda^{1/d}}{d}t^{(1-\nu d}$,and, for the case of the size-dependent migration process, it is

(14)

143

$\overline{V_{t}}\propto\frac{\mu^{1/d}}{d}\cdot N(t)$.

$[ \int_{0^{t}}N(\tau)d\tau]^{1-dy_{d}}$

.

(15)

LOGISTIC GROWTH OF POPULATIONAs an example of the size-dependent migration process, we deal with the casewhen the total population size of the group of colonies grows in the logisticmanner (Fig. $3(a)$ ):$N(t)=N(0) \cdot\{(1-\frac{N(0)}{K})e^{-}’+\frac{N(0)}{K}\backslash /^{-1}$

(16)

,

where is the intninsic growth rate of the population and$\epsilon$$K$

is the carrying

capacity for the total population of the group. In this case, some fundamental
calculations show
(17)$k)_{l}= \frac{\mu}{\epsilon}\cdot K\cdot\ln\{\frac{N(0)}{K}\cdot(\not\in-1)+1\}$$\#k=\frac{1}{\epsilon}\ln\{\frac{K}{N(0)}\cdot(e^{(\epsilon 1\downarrow\downarrow)k/\kappa_{-1)+}}1\}.$(18)

Generic feature of these expected values is shown in Fig. 3(b) and Fig. 3(c). The
time interval between the nearest two colony productions is given by$\#k+1-\#k=_{\frac{1}{\epsilon}}\ln\{\frac{e^{\epsilon/(\mathscr{O}}+\frac{N(0)}{K}(1-\frac{N(0)}{K})\cdot e^{-(\epsilon 1\downarrow\downarrow)k/K}}{1+\frac{N(0)}{K}(1-\frac{N(0)}{K})\cdot e^{-(\epsilon 1\downarrow\iota)k/K}}\{$.
This value tends to a constant $l/pK$ as$karrow\infty$(19), which means that the colony

production is expected to occur periodically. In addition, from (17) and (18), forsufficiently large and sufficiently large ,
$t$

$k$

144

$N\cup\dagger 1$

BER OF COLONIES

Fig. 3. The size-dependent immmigration model for the logistic growing population.
(a) Typica1 time-development of logistic population growth; (b) Typical $\#\#:(c)$Typicallk$\cdot$(20)$\phi h^{\sim}t$$\#k^{\sim k}$

.

(21)

That is, these expected values increase linearly in a sufficiently grown group,
which is the same as for the case of the free migration process.
In this case, the expanding way of settlement range essentially depends onthe fractal dimension (Fig. $4(a,$ ). The expected velocity of expanding ofsettlement range is expressed as follows:
$d$

$b)$

$\overline{V_{t}}\propto\frac{\epsilon}{d}(\frac{\mu K}{\epsilon})^{1/d}\frac{N(0)}{K}\cdot\frac{e^{\sigma}}{\{\frac{N(0)}{K}\cdot(\theta-1)_{/}+1^{\backslash }[\ln\{\frac{N(0)}{K}\cdot(\theta-1)+1^{1_{(}}]^{1-11d}}$

For $1<d\leq 2$ , this expected velocity decreases to zero at a sufficiently large time
(Fig. $4(c,$ ). This means that, for a sufficiently grown group, the velocity of the
expanding of settlement range is very small, while the number of colonies
continuously increases; that is, the new colonies tend to be produced within thevacant areas among the pre-settled colonies. On the other hand, the timedevelopment of the expected velocity in the earlier period depends on the initial
population size of the group (Fig. $4(d)$ ). In the case when the initial population is
sufficiently large, the expected velocity monotonically decreases in time, while inthe case when it is small, the expected velocity increases in the earlier period anddecreases after peaking. Analytically, if the following condition is satisfied, theformer case occurs, and otherwise the latter (Fig. $4(d)$ ):$d)$$1-(1- \frac{1}{d})\cdot e^{1/d}\leq\frac{N(0)}{K}$

.

(23)

.

(22)

145

Fig. 4. In the size-dependent immigration process model for the logistic population
growth, the contnibution of the fmctal dimension and the initial population size
$N(O)/K$ to: $(a, b)$ the time-development of the expected settlement range , that is,
$d$

$R$

$kt^{/d}$: $(c, d)$ the time-development of the expected velocity 17. For (c), $N(O)/K=$

0.1, and for (d), $d=1.5$ . The graph shape of the expected velocity
and $N(O)/K$ as shown totally in the figure attached to (d).

$\overline{V_{t}}$depends on$d$This condition can be easily derived by examining the sign of the t-derivative of(22).

In the period when the population size $N(t)$ is sufficiently small, thepopulation growth can be well-approximated by exponential growth. In thisperiod, the same argument as above gives the approximate results on the behavior

of colony dispersal as follows:$N(t)\approx N(0ffl$(24)$k h^{\approx}\frac{\mu}{\epsilon}N(0X\not\in-1)$(25)$\#k^{\approx}\frac{1}{\epsilon}\ln(\frac{\epsilon/\mu}{N(0)}k+1)$(26)$\#k+1-\#k^{\approx}\frac{1}{\epsilon}\ln(\frac{N(0)+\epsilon(k+1y\mu}{N(0)+\epsilon H\mu})$(27)$\overline{V_{t}}\propto\frac{\epsilon}{d}(\frac{\mu N(0)}{\epsilon})^{1/d}\frac{e^{\epsilon r}}{(e^{\sigma}-1)^{1-11d}}$

.

(28)

146

屋
$\underline{\vee\vee-\varpi=}$longitude

Fig. 5. 106 neolithic European sites in the geographic coordinates, used by
Ammerman and Cavalli-Sforza (1984). The black square indicates the oldest site,
Aswad (9690 B.C., C-14 date). Black discs are for those sites before 5800 B.C. (C14 date), and white ones for those after 5800 B.C.

In this case, the expected velocity increases exponentially while the number ofcolonies grows exponentially.

SPREAD OF EARLY FARMING IN EUROPE
Ammerman and Cavalli-Sforza (1984) calculated the isochron map of the spreadof early farming in Europe from the data of 106 neolithic European sites with C14 dates (9690 B.C. -4160 B.C.). The computer-generated isochron map gives
the impression that early farming might have spread in a spatially continuous
manner in Europe. This is an approximation to the spatial spread through theanalogy of diffusion process. However, in contrast to the spatial spread of various
species of animals, insects, and plants, the spatial spread of a group of humans
frequently involves the production spatially disconnected units, that is, colonies.
The spatial distribution expands essentially by a series of productions of newcolonies. The spread of early farming in Europe, dealt with by Ammerman andCavalli-Sforza (1984), can be regarded as such a case.
In this section, we apply our mathematical model described above to thedata and estimate the parameters of the model to try to discuss some features ofthe spread of early farming in Europe.
As for the way of population growth, we assume the exponential one givenby (24). This is appropriate in the case when the population growth does notcause the depletion of environmental factors necessary for survival.

147

time

Fig. 6. Time-development of the number of neolithic sites, which is cumulated after
the oldest site, Aswad (9690 B.C., C-14 date). Time axis shows the C-14 datepassed $aRer$ 9690 B.C. All 106 neolithic sites are plotted for the data of Ammerman
and Cavali-Sforza (1984). $k1$ curve for the exponential growth in the sizedependent migration process model is overlaid, fit to 76 data of neolithic sites$($ before 5800 B.C. (C-14 date). $\mu N(O)/\epsilon=3.352$ and $\epsilon=8.233x10^{-4}$ in(25).

Fig. 5 shows the 106 neolithic European sites in the geographic
coordinates used by Ammerman and Cavalli-Sforza (1984). Their spatial
distribution seems to show inhomogeneity in direction. Beginning with the oldest
site, Aswad (9690 B.C.; 33. $36N,$ $36.30E$), we count the cumulative number ofcolonies in order of descending C-14 date as shown in Fig. 6. Plots in the figure
indicate that the continuity of the time-development of the number of colonies
seems to break at around 3800 years after Aswad (i.e., around 5800 B.C.). Thus,
we use only the 76 data before 5800 B.C., up to the site Reichtett (5940 B.C.;

48. $6N,$ $7.75E$).
Since Fig. 6 can be regarded as coIresponding to the time-development of, we try to fit } given by (25) to the data. The result is overlaid in Fig. 6. The
estimated parameters result in $\mu N(O)/\epsilon=3.572;\epsilon=8.233^{x}10^{-4}$ .

$k)_{t}$$\ell_{\{}$$i$

Next, we try to estimate the ffactal dimensionthat characterizes theof the settlement area and$d$pattern of spatial distribution. From (11), the range
the number $M$ of colonies within it have the relation: $\log M=d\cdot\log R+const$.
Therefore, we can estimate from the slope of the line fit to the plots of log $M$
against $\log R$ . The diameter can be calculated from the data of the locations ofneolithic sites (Fig. 7). We use the gyration-radius method to estimate theparameter (as for the method, see, for instance, Mandelbrot, 1982). All 76 sites
before 5800 B.C. are considered. The number of sites distributed within the disccentered at the oldest site, Aswad, is counted. For disc radius large enough tocontain more than 10 colonies, the plot of the number of sites against the disc
$R$

$d$

$d$

148

time

Fig. 7. Expanded diameter of the settlement range. Time axis shows the C-14 datepassed after 9690 B.C.

$\underline{\frac{\wedge oc\wedge}{*\vee oooo}}$log(r)

Iog(expanded diameter)

Fig. 8. (a) Number of sites within the distance from the oldest site, Aswad, in loglog coordinates. For the distance that contains more than 10 sites, plots fit well tothe line with the slope 1.671. (b) Relation between the expanded diameter ofsettlement area and the number of colonies in log-log coordinates. The overlaid line
indicates the slope 1.671. The unit of measured distance is conventionally selected.
$r$

$r$

radius can be fitted well by a straight line with slope 1.671, as estimated by theleast-square method (Fig. $8(a)$ ). Hence, the spatial distribution of neolithic sites is
estimated to have the characteristic fractal dimension $d=1.671$ . Since thediameter and the number of neolithic sites are time-dependent, it is likely that thepaIameter might change in time in the period considered now. However, as Fig.
8(b) shows, the time-dependent relation between the expanded diameter and thenumber of colonies in log-log coordinates, the estimated $d=1.671$ even holds
well. Therefore, we deal with as time-independent constant: $d=1.671$ .
From (28) with these estimated parameters, the time-development of the
expected velocity of the expanding of settlement range can be drawn, resulting in
Fig. 9(a). It can be seen that the velocity is relatively small in the first century
after Aswad and then increases exponentially. In Fig. 9(b), the same expected
$d$

$d$

149

time

#colony

Fig. 9. Time-development of the expected velocity (28) of the expanding of thesettlement area, for $\mu N(O)/\epsilon=3.352;\epsilon=8.233x10^{-4};d=1.671$ . $(a)$ timedevelopment of the expected velocity (28); (b) the expected velocity (28) against the
expected number of colonies (25).

velocity is plotted against the time-development of the expected number ofcolonies, which follows (25). As the number of colonies becomes sufficiently
large, the velocity of the expanding of the settlement area increases.
REFERENCES
Ammerman, A. J. and Cavalli-Sforza, L. L. (1984) Neolithic Transition and The Genetics ofPopulations in Europe, Princeton University Press, Princeton, New Jersey.
Bartlett, M. S., F.R.S. (1978) An Introduction to Stochastic Processes, Cambridge University
Press, Cambridge.
Britton, N. F. (1986) Reaction-Diffision Equations and Their Applications to Biology, Academic
Press, London.
Mandelbrot, B. B. (1982)
Fractal Geometry ofNature, Freeman, San Francisco.
Martin, P. S. (1973) The discovery of America. Science. 179, 969-974.
Mosinam, J. E. and Martin, P. S. (1975) Simulating overkill by paleoindians. Am SCL 63, 304$Th\ell$313.

Murray, J. D. (1989) Mathematical Biology, Springer-Verlag, New York.
Okubo, A. (1980) Diffision and Ecological Problems: Mathematical Models, Springer-Verlag,
New York.
Skellam, J. G. (1951) Random dispersal in theoretical populations. Biometrika. 38, 196-218.