! .../Projects/2025_Perez_real_numbers/coast_Fort/37a_de_28_Fig_4__eta_0/2026_xx_xx__Fig_3.f90 ! Author: J. Fort ! This code is a modification of the code Fort, AAS (2022). The major change is coupled logistic ! reproduction here (vs. indep. logistics in AAS 2022). Genetics as in Sci.Rep.2017 & Nat.Com.2025. !---------------------------------------------------------------------------------------- integer::i,j,y,t,z,x,n,m,Cyprus,CAnatolia,DAnatolia,EAnatolia,Crete,GGreece real,dimension (6999,6999)::k,h,a,an,ap,bn,bp real,dimension (200)::rfront real,dimension (6999)::East,Eastm,previous,ee real dis,pe,rr,am,rp,bm,f,ga,fraction ! Iterations = final number of generations: t=172 because Data S2, region 5, 171 gen. t=172 ! Horizontal (or vertical, for ga=1 and eta=f) cultural transmission parameters: f=0 ga=1 ! distance between nodes (km): dis=1 ! Ro=2.45 as in Fort, PNAS 2017. In fortran LOG is the natural logarithm, so ! a=rr=ln(Ro)=ln(2.45)=0.896 gen^-1 ˇ 1 gen/32 yr =0.028 yr^-1: rr=LOG(2.45) ! Mesolithic initial growth rate in gen^-1 (Fort, Pujol & Cavalli-Sforza 2004): rp=0.59 ! Neolithic satutarion density (Currat+Excoffier 2005): am=1.28 ! Mesolithic saturation density (Currat+Excoffier 2005): bm=0.064 ! Persistency (Fort et al., Phys. Rev. E 2007): pe=0.38 !================================================================================ ! INITIAL CONDITIONS: fraction= 0.474 ! Fraction of farmers with haplogroup K at point O. ! Initally only farmers at point O (y=100); an=K farmers, ap=non-K farmers, k=HGs: an=0 ap=0 a=0 an(1,100)=1.28*fraction ap(1,100)=1.28-an(1,100) a(1,100)=an(1,100)+ap(1,100) bn=0 bp=0 ! Initially HGs everywhere except at point O: k=0.064 k(1,100)=0 h=0 ! We keep the initial Neolithic profile: OPEN (5,FILE='NEOPROF0.DAT',STATUS='UNKNOWN') do zz=1,6999 WRITE(5,*) zz,a(1,zz) end do CLOSE(5) ! We keep the initial Mesolithic profile: OPEN (5,FILE='MESOPROF0.DAT',STATUS='UNKNOWN') do zz=1,6999 WRITE(5,*) zz,k(1,zz) end do CLOSE(5) ! This deletes the file RFRONT.DAT: OPEN (2,FILE='RFRONT.DAT',STATUS='UNKNOWN') WRITE (2,*) ' ' CLOSE (2) !================================================================================ ! 1) INTEGRATION IN TIME: do j=1,t !================================================================================ ! The dispersal in the WESTERN Mediterranean begins as close as possible to the intersection of the 2 regresssions (121.5 gen and 1411 km from point O). ! We know from running this code that at t=121 gen rfront(121)=1394 km, i.e. the EASTERN front is at 1394 km from point O or y=1394+100=1494 km. ! We use the %K from the EASTERN expansion as initial %K for the WESTERN expansion, but there are no farmers at y=1494 km because 1394/18=77.44 is not ! an integer. The closest node to y=1494 km (1394 km from point O) to which farmers arrive is y=1504 km (1404 km from point O) because 1404/18=78 is ! an integer. So we have to evaluate the %K at y=1504 km and apply it at y=1494 km, immediately after all calculations for the EASTERN expansion and ! t=121 gen, so t=122 gen: if (j==122) then !print *, an(1,1504) ! result: 1.990595E-02, different from zero-->OK !print *, ap(1,1504) ! result: 2.212932E-02, different from zero-->OK bcd=an(1,1504) efg=ap(1,1504) an=0 ap=0 a=0 an(1,1494)=bcd*1.28/(bcd+efg) ! This line and next imply an+ap=1.28=am, and this is necessary to obtain rfront.dat and Fig1b.dat as for Fig. 1b. ap(1,1494)=efg*1.28/(bcd+efg) bn=0 bp=0 b=0 k=0.064 k(1,1494)=0 h=0 endif ! 2) REPRODUCTION: COUPLED LOGISTIC TERMS FROM LAPOLICE, WILLIAMS & HUBER, Nature Comm. 2025: do z=1,6999 do x=1,6999 ! We neglect nodes with approx. zero population density (otherwise it is very slow): if (an(x,z)<1D-10) then else bn(x,z) = an(x,z) + rr * an(x,z) * (1- ((an(x,z)+ap(x,z)+k(x,z))/am) ) endif if (ap(x,z)<1D-10) then else bp(x,z) = ap(x,z) + rr * ap(x,z) * (1- ((an(x,z)+ap(x,z)+k(x,z))/am) ) endif if (k(x,z)<1D-10) then else k(x,z) = k(x,z) + rp * k(x,z) * (1- ((an(x,z)+ap(x,z)+k(x,z))/bm) ) end if an(x,z) = bn(x,z) ap(x,z) = bp(x,z) end do end do bn=0 bp=0 !================================================================================= ! 3) INTERBREEDING: h=k do z=1,6999 do x=1,6999 ! We neglect nodes with approx. zero Mesol. popul. density (otherwise, floating-point error!): if (k(x,z)<1E-10) then else k(x,z)=h(x,z)-f*h(x,z)*(an(x,z)+ap(x,z))/(an(x,z)+ap(x,z)+ga*h(x,z)) ap(x,z)=ap(x,z)+f*h(x,z)*(an(x,z)+ap(x,z))/(an(x,z)+ap(x,z)+ga*h(x,z)) end if if (k(x,z)<0) then k(x,z)=0 end if if (a(x,z)>1.28) then a(x,z)=1.28 end if end do end do h=0 !================================================================================ ! 4) NEOLITHIC DISPERSAL: ! The dispersal in the WESTERN Mediterranean begins at 122 generations, see also line 200 below: if (j>121) then goto 200 endif ! ======DISPERSAL IN THE EASTERN MEDITERRANEAN:================================== ! length of coastal jumps (km) in the EASTERN Mediterranean: n=18 ! length of inland jumps (km) in the EASTERN Mediterranean: m=18 ! INLAND jumps: i=2,6999 (not i=1,6999) because the coast is the straight line (i,y) with i=1: do y=1,6999 do i=2,6999 !----------------------- ! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!): if (an(i,y)+ap(i,y)<1E-10) then goto 100 endif !----------------------- do z=-1,1 do x=-1,1 ! a(i,y) is the density of farmers at the initial node of jumps: if (x==0.and.z==0) then bn(i,y)=bn(i,y)+(pe)*an(i,y) bp(i,y)=bp(i,y)+(pe)*ap(i,y) goto 80 endif if(x==0.or.z==0) then bn(i+m*x,y+m*z)=bn(i+m*x,y+m*z)+((1-pe)/4)*an(i,y) bp(i+m*x,y+m*z)=bp(i+m*x,y+m*z)+((1-pe)/4)*ap(i,y) goto 80 endif 80 continue end do end do 100 end do end do !*********************** ! COASTAL jumps: the coast is the straight line (i,y) with i=1: do y=1,6999 !----------------------- ! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!): if (an(1,y)+ap(1,y)<1E-10) then goto 150 endif !----------------------- do z=-1,1 do x=-1,1 if (x==0.and.z==0) then bn(1,y)=bn(1,y)+(pe)*an(1,y) bp(1,y)=bp(1,y)+(pe)*ap(1,y) goto 120 endif ! From coastal nodes, farmers jump to 3 nodes (not to 4 as above line 100): if(x==0.and.z==-1) then bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y) bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y) goto 120 endif if(x==0.and.z==1) then bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y) bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y) goto 120 endif ! this is an inland jump, so its lenght is m, not n: if(x==1.and.z==0) then bn(1+m*x,y)=bn(1+m*x,y)+((1-pe)/3)*an(1,y) bp(1+m*x,y)=bp(1+m*x,y)+((1-pe)/3)*ap(1,y) goto 120 endif ! Node (x=-1.and.z==0) is not included because (1+x,y+z)=(0,y) is sea (not coast). 120 continue end do end do 150 end do goto 370 ! ======DISPERSAL IN THE WESTERN MEDITERRANEAN:============================== ! length of coastal jumps (km) in the WESTERN Mediterranean: 200 n=430 ! length of inland jumps (km) in the WESTERN Mediterranean: m=50 ! The dispersal in the WESTERN Mediterranean begins as close as possible to the intersection of ! the 2 regresssions. We know from running this code that at t=121 gen the EASTERN front is at ! rfront=1394 km from point O, i.e. y=1394+100=1494 km. So we assume that the WESTERN front starts ! at this position and 122 gen: ! INLAND jumps: i=2,6999 (not i=1,6999) because the coast is the straight line (i,y) with i=1: do y=1494,6999 do i=2,6999 !----------------------- ! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!): if (an(i,y)+ap(i,y)<1E-10) then goto 300 endif !----------------------- do z=-1,1 do x=-1,1 ! a(i,y) is the density of farmers at the initial node of jumps: if (x==0.and.z==0) then bn(i,y)=bn(i,y)+(pe)*an(i,y) bp(i,y)=bp(i,y)+(pe)*ap(i,y) goto 280 endif if(x==0.or.z==0) then bn(i+m*x,y+m*z)=bn(i+m*x,y+m*z)+((1-pe)/4)*an(i,y) bp(i+m*x,y+m*z)=bp(i+m*x,y+m*z)+((1-pe)/4)*ap(i,y) goto 280 endif 280 continue end do end do 300 end do end do !*********************** ! COASTAL jumps: the coast is the straight line (i,y) with i=1: do y=1494,6999 !----------------------- ! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!): if (an(1,y)+ap(1,y)<1E-10) then goto 350 endif !----------------------- do z=-1,1 do x=-1,1 if (x==0.and.z==0) then bn(1,y)=bn(1,y)+(pe)*an(1,y) bp(1,y)=bp(1,y)+(pe)*ap(1,y) goto 320 endif ! From coastal nodes, farmers jump to 3 nodes (not to 4 as above line 100): if(x==0.and.z==-1) then bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y) bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y) goto 320 endif if(x==0.and.z==1) then bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y) bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y) goto 320 endif ! this is an inland jump, so its lenght is m, not n: if(x==1.and.z==0) then bn(1+m*x,y)=bn(1+m*x,y)+((1-pe)/3)*an(1,y) bp(1+m*x,y)=bp(1+m*x,y)+((1-pe)/3)*ap(1,y) goto 320 endif ! Node (x=-1.and.z==0) is not included because (1+x,y+z)=(0,y) is sea (not coast). 320 continue end do end do 350 end do 370 continue !********************** ! We keep the Neolithic population in matrices a, to be used in the next iteration: an=bn ap=bp ! We reinitiate matrices b, they will change during the next iteration: bn=0 bp=0 a=an+ap !=============================================================================== ! 5) WE KEEP SOME RESULTS BEFORE JUMPING TO THE NEXT TIME STEP: !------------------ ! 5a) FRONT POSITION along the vertical direction as a function of time: ! the slope is the spread rate in km/generation (import with e.g. origin + linear fit) OPEN (2,ACCESS='APPEND',FILE='RFRONT.DAT',STATUS='UNKNOWN') ! I idenfity the first node with rho<0.1, and find the point with rho=0.1 (by fitting a straight line): ! COAST=EDGE OF THE GRID =(1,zz) [I use zz because z causes error "cannot link"!]: do zz=6999,100,-1 ! the upper limit of coast (x,y)=(1,6999) fails, so I look at y<6250: if (zz>6250) then goto 400 endif ! zz+n appears below because the node zz+1 is always empty of farmers: they arrive only to coastal nodes separated n km. if (a(1,zz)>0.1) then rfront(j)=dis*(zz+(a(1,zz)-0.1)*n/(a(1,zz)-a(1,zz+n))-100) exit end if 400 continue end do IF (.not.(j<1)) WRITE (2,*) j,rfront(j) close(2) !------------------ ! 5b) FRONT PROFILES: ! The dispersal in the WESTERN Mediterranean begins as close as possible to the intersection of ! the 2 regressions. We know from running this code that at t=121 gen the EASTERN front is at y=1394+100=1494 km ! Neolithic profile at generation 140: if (j==140) then OPEN (4,FILE='NEOPROF.DAT',STATUS='UNKNOWN',POSITION='APPEND') do zz=1,1493 WRITE(4,*) zz,East(zz) end do do zz=1494,6999 WRITE(4,*) zz,a(1,zz) end do CLOSE(4) end if ! Mesolithic profile at generation 140: if (j==140) then OPEN (5,FILE='MESOPROF.DAT',STATUS='UNKNOWN',POSITION='APPEND') do zz=1,1493 WRITE(5,*) zz,Eastm(zz) end do do zz=1494,6999 WRITE(5,*) zz,k(1,zz) end do CLOSE(5) end if !------------------ ! 5c) Dates for Fig. 1b (using distances along the coast in km from Data S1): OPEN (6,FILE='FIG1b.DAT',STATUS='UNKNOWN',POSITION='APPEND') ! Eastern Mediterranean: if (j==1) then WRITE(6,*) 0,0 endif ! Below, no need to subtract 100 km (initial position of farmers) because already subtracted in rfront above. if (Cyprus==1) then goto 500 endif if (rfront(j)>453) then WRITE(6,*) 453,j-1+(453-rfront(j-1))/(rfront(j)-rfront(j-1)) Cyprus=1 endif 500 continue if (CAnatolia==1) then goto 503 endif if ((rfront(j))>583) then WRITE(6,*) 583,j-1+(583-rfront(j-1))/(rfront(j)-rfront(j-1)) CAnatolia=1 endif 503 continue if (DAnatolia==1) then goto 510 endif if ((rfront(j))>846) then WRITE(6,*) 846,j-1+(846-rfront(j-1))/(rfront(j)-rfront(j-1)) DAnatolia=1 endif 510 continue if (EAnatolia==1) then goto 513 endif if ((rfront(j))>1122) then WRITE(6,*) 1122,j-1+(1122-rfront(j-1))/(rfront(j)-rfront(j-1)) EAnatolia=1 endif 513 continue if (Crete==1) then goto 516 endif if ((rfront(j))>1172) then WRITE(6,*) 1172,j-1+(1172-rfront(j-1))/(rfront(j)-rfront(j-1)) Crete=1 endif 516 continue if (GGreece==1) then goto 520 endif if ((rfront(j))>1231) then WRITE(6,*) 1231,j-1+(1231-rfront(j-1))/(rfront(j)-rfront(j-1)) GGreece=1 endif 520 continue if (j==121) then WRITE(6,*) NINT(rfront(121)),121 endif ! Western Mediterranean: ! Using the distance of the oldest site per region does not lead to a straight line. I use one node every n km: 690 continue do xx=4,14 if (ee(xx)==1) then goto 700 endif if (a(1,1494+n*xx)>0.1) then WRITE(6,*) 1494+n*xx-100,j-(a(1,1494+n*xx)-0.1)/(a(1,1494+n*xx)-previous(1494+n*xx)) ee(xx)=1 end if 700 continue end do CLOSE(6) !=============================================================================== ! This is used to obtain the times in Fig. 1b (otherwise we get no straight line in the West Med.): do yy=1,6999 previous(yy)=a(1,yy) end do !------------------ !5d) Numbers of generations, km from point O and %K for Fig. 3 (regions in Data S2): OPEN (7,FILE='FIG3.DAT',STATUS='UNKNOWN',POSITION='APPEND') !--------------------------------------------------------------------- ! We include some values of %K in Fig.3 between point 0 and region 1 because it is a curve, not a straight line. ! The time (j) is not important because %K changes <1% in 20 gen (Fig. S6c in Fort & Pérez, Nature Comm. 2025). ! We need a node with farmers: its distance from point O is an integer multiple of n. y is this distance +100km. if (j==50) then y = 17*n+100 write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y)) endif if (j==70) then y = 34*n+100 write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y)) endif if (j==90) then y = 50*n+100 write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y)) endif if (j==120) then y = 67*n+100 write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y)) endif if (j==121) then y = 77*n+100 write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y)) endif !--------------------------------------------------------------------- selectcase(j) case(1) ! 1 is the number of generations j; below, j=0 gen and i=0 km for region 1 = initial region (it contains point O): write(7,*) 0, 0, 47.4 !--------------------------------------------------------------------- case(126) ! Region 2 = Greece and North Macedonia: 126 generations. In fact Data S2 says 116 generations, so in ! principle this should be case(116) = 116 generations. However, then "Error M6101. Flotaing-point error" ! because the front has not yet arrived at 1670+100 km (Fig 1b or rfront.dat), i.e. an=0, ap=0 so %K=0/0. ! It is reasonable that it has not arrived because the model is homogeneous. ! But we know that at 126 gen it has arrived from rfront.dat. Using a larger value of time ! is OK because %K changes <1% in 20 generations (Fig. S6c in Fort & Pérez-Losada, Nature Comm. 2025). y = 1670+100 ! 100 km is the initial position of farmers, i.e. the position of point O in this code. ! 1670 km is the distance from point O to region 2 (Data S2). ! The node y=1670+100 km is empty of farmers because they jump 430 km from y=1494 km, so an=0, ap=0 so %K=0/0. ! We look for the closest node where farmers have arrived from a coastal jump: abc=ABS(y-1494) do xx=1,14 if (ABS(y-(1494+n*xx))