Lizard evolution simulations

Setup

<< Statistics`DiscreteDistributions`

    
<< Graphics`Graphics`

    
data  =  
Table[
Random[
            PoissonDistribution[10]], 
{i, 1, 
            2000}
];
\(\(Histogram[
data, 
AxesStyle\  -> \ {Thickness[
              0.005]}, 
TextStyle -> {FontWeight -> "Bold", 
            FontSize -> 
              18}
];\)
\)

Spatial distribution of lizards with varying types

Mo0  =  25;
Molist  = 
Table[
{Random[], 
            Random[]}, 
{i, 1, 
            Mo0}
];
Mb0  =  50;
Mblist  = 
Table[
{Random[], 
            Random[]}, 
{i, 1, 
            Mb0}
];
My0  =  25;
Mylist  = 
Table[
{Random[], 
            Random[]}, 
{i, 1, 
            My0}
];
F0  =  200;
Flist  = 
Table[
{Random[], 
            Random[]}, 
{i, 1, 
            F0}
];

    
\(\(\(draw[lizards_, 
          color_]\  := 
Module[{r, g, 
            b}, 
{r, g, b}\  = \ 
            color; 
If[
            Length[lizards] > 
              0, 
Graphics[
{RGBColor[r, 
                  g, b], 
Table[
Point[
                    lizards[\([i]\)]], 
{i, 1, 
                    Length[lizards]}
]
}
], 
Graphics[
{RGBColor[1, 1, 
                  1], 
PointSize[
                  0], 
Point[{0, 
                    0}]
}
]
]
];\)\(
\)
    \)\)
Show[
draw[
          Molist, {1, 0.25, 0}], 
draw[
          Mblist, {0, 0, 1}], 
draw[
          Mylist, {1, 1, 0}], 
draw[
          Flist, {0, 0, 0}], 
Prolog -> 
          PointSize[
            0.015], 
PlotRange -> {{-0.01, 
              1.01}, {-0.01, 1.01}}, 
Axes -> 
          False, 
AspectRatio -> 
          Automatic, 
DisplayFunction -> $DisplayFunction
];

Moving lizards

distance[{x1_, y1_}, {x2_, y2_}]  := 
Module[{x, 
            y}, 
x  =  
            Min[{Abs[x1 - x2], Abs[x1 - x2 + 1], 
                Abs[x1 - x2 - 1]}]; 
y  =  
            Min[{Abs[y1 - y2], Abs[y1 - y2 + 1], 
                Abs[y1 - y2 - 1]}]; 
Sqrt[
            x^2  +  y^2]
];

    
ro  =  0.2;
rb  =  ro/2;

    
conflict[Molist_, Mblist_, {x_, y_}, 
          rp_]  := 
Module[{i}, 
Sum[
E^(-(distance[
                          Molist[[i]], {x, y}]/(ro + 
                            rp))), 
{i, 1, 
                Length[Molist]}
]  + 
Sum[
E^(-(distance[
                          Mblist[[i]], {x, y}]/(rb + 
                            rp))), 
{i, 1, 
                Length[Mblist]}
]
];

    
δ  =  0.01;
step[{x_, 
            y_}] := 
Module[{θ}, 
θ  =  
            2*Pi*Random[]; 
{Mod[
              x + δ* Cos[θ], 1], 
Mod[
              y + δ* Sin[θ], 
              1]
}
];
bigδ  =  0.1;
bigstep[{x_, 
            y_}]  := 
Module[{θ}, 
θ  =  
            2*Pi*Random[]; 
{Mod[
              x + bigδ* Cos[θ], 1], 
Mod[
              y + bigδ* Sin[θ], 
              1]
}
];

    
pdeath  =  1/26;
pbdeath  =  1/10;

    
move[Molist_, Mblist_, Mylist_, Flist_, babyMolist_, babyMblist_, 
          babyMylist_, 
          babyFlist_]  := 
Module[
{i, 
            newMolist, newMblist, newMylist, 
            newFlist, 
newbabyMolist, newbabyMblist, 
            newbabyMylist, 
            newbabyFlist
}, 
newMolist  = {}; 
Do[
If[
              Random[]  >  pdeath, 
newposition  =  
                step[Molist[[i]]]; 
If[
                conflict[Delete[Molist, i], Mblist, newposition, 
                    ro]  < 
conflict[Delete[Molist, i], 
                    Mblist, Molist[[i]], ro], 
AppendTo[
                  newMolist, newposition], 
AppendTo[
                  newMolist, 
                  Molist[[i]]]
]
], 
{i, 1, 
              Length[Molist]}
]; 
newMblist  =  {}; 
Do[
If[
              Random[]  >  pdeath, 
newposition  =  
                step[Mblist[[i]]]; 
If[
                conflict[Molist, Delete[Mblist, i], newposition, 
                    rb]  < 
conflict[Molist, 
                    Delete[Mblist, i], Mblist[[i]], 
                    rb], 
AppendTo[newMblist, 
                  newposition], 
AppendTo[newMblist, 
                  Mblist[[i]]]
]
], 
{i, 1, 
              Length[Mblist]}
]; 
newMylist  =  {}; 
Do[
If[
              Random[]  >  pdeath, 
newposition  =  
                step[Mylist[[i]]]; 
AppendTo[newMylist, 
                newposition]
], 
{i, 1, 
              Length[Mylist]}
]; 
newFlist  =  {}; 
Do[
If[
              Random[]  >  pdeath, 
newposition  =  
                step[Flist[[i]]]; 
AppendTo[newFlist, 
                newposition]
], 
{i, 1, 
              Length[Flist]}
]; 
newbabyMolist  =  {}; 
Do[
If[
              Random[]  >  pbdeath, 
newposition  =  
                bigstep[babyMolist[[i]]]; 
AppendTo[
                newbabyMolist, 
                newposition]
], 
{i, 1, 
              Length[babyMolist]}
]; 
newbabyMblist  =  {}; 
Do[
If[
              Random[]  >  pbdeath, 
newposition  =  
                bigstep[babyMblist[[i]]]; 
AppendTo[
                newbabyMblist, 
                newposition]
], 
{i, 1, 
              Length[babyMblist]}
]; 
newbabyMylist  =  {}; 
Do[
If[
              Random[]  >  pbdeath, 
newposition  =  
                bigstep[babyMylist[[i]]]; 
AppendTo[
                newbabyMylist, 
                newposition]
], 
{i, 1, 
              Length[babyMylist]}
]; 
newbabyFlist  =  {}; 
Do[
If[
              Random[]  >  pbdeath, 
newposition  =  
                bigstep[babyFlist[[i]]]; 
AppendTo[
                newbabyFlist, 
                newposition]
], 
{i, 1, 
              Length[babyFlist]}
]; 
Show[
draw[
              newMolist, {1, 0.25, 0}], 
draw[
              newMblist, {0, 0, 1}], 
draw[
              newMylist, {1, 1, 0}], 
draw[
              newFlist, {0, 0, 0}], 
draw[
              newbabyMolist, {0, 1, 0}], 
draw[
              newbabyMblist, {0, 1, 0}], 
draw[
              newbabyMylist, {0, 1, 0}], 
draw[
              newbabyFlist, {0, 1, 0}], 
Prolog -> 
              PointSize[
                0.015], 
PlotRange -> {{-0.01, 
                  1.01}, {-0.01, 1.01}}, 
Axes -> 
              False, 
AspectRatio -> 
              Automatic, 
DisplayFunction -> $DisplayFunction
]; 
 (*
all  =  
              Union[newMolist, 
                newMblist]; 
ave  =  
Sum[
distance[all[[i]], 
                    all[[j]]], 
{i, 1, 
                    Length[all] - 1}, {j, i + 1, 
                    Length[all]}
]/(Length[
                      all]*(Length[all] - 1)); 
Print[
              ave];
*) 
{newMolist, 
            newMblist, newMylist, newFlist, 
newbabyMolist, 
            newbabyMblist, newbabyMylist, 
            newbabyFlist
}
];


Do[
{Molist, Mblist, Mylist, 
            Flist, 
babyMolist, babyMblist, babyMylist, 
            babyFlist
}  =  
move[
Molist, Mblist, Mylist, 
            Flist, 
babyMolist, babyMblist, babyMylist, 
            babyFlist
], 
{t, 1, 
          26}
];

Mating of lizards

pmate  =  0.25;
aveobabies  =  10;
avebbabies  =  10;
aveybabies  =  10;
ry  =  (4/3)*
          rb;

mate[Molist_, Mblist_, Mylist_, 
          Flist_]  := 
Module[
{babyMolist, babyMblist, 
            babyMylist, babyFlist, 
Nobabies, Nbbabies, 
            Nybabies, 
bmates, bstops, ymates, i, 
            j}, 
babyMolist  =  {}; 
babyMblist  =  {}; 
babyMylist  =  {}; 
babyFlist  =  {}; 
Do[
If[
              Random[] < 
                pmate, 
omates  =  {}; 
bmates  =  {}; 
ymates  =  {}; 
ystoppers  =  {}; 
Do[
dj  = 
                  distance[Flist[[i]], 
                    Molist[[j]]]; 
If[
                  dj < ro, 
AppendTo[
                    omates, {j, 
                      dj}]
], 
{j, 1, 
                  Length[Molist]}
]; 
Do[
dj  = 
                  distance[Flist[[i]], 
                    Mblist[[j]]]; 
If[
                  dj < rb, 
AppendTo[
                    bmates, {j, 
                      dj}]
], 
{j, 1, 
                  Length[Mblist]}
]; 
Do[
dj  = 
                  distance[Flist[[i]], 
                    Mylist[[j]]]; 
If[
                  dj < ry, 
AppendTo[
                    ymates, {j, 
                      dj}]
], 
{j, 1, 
                  Length[Mylist]}
]; 
If[
                Length[omates] ==  1  &&  
                  Length[bmates] ==  0  &&  
                  Length[ymates] ==  
                    0, 
Nobabies  =  
                  Random[PoissonDistribution[
                      
                      aveobabies]]; 
Do[
If[
                    Random[] < 1/2, 
AppendTo[babyMolist, 
                      Flist[[i]]], 
AppendTo[babyFlist, 
                      Flist[[i]]]
], 
{b, 1, 
                    Nobabies}
]
]; 
If[
                Length[bmates] ==  1  &&  
                  Length[omates] ==  
                    0, 
Nbbabies  =  
                  Random[PoissonDistribution[
                      avebbabies]]; 
Do[
If[
                    Random[] < 1/2, 
AppendTo[babyMblist, 
                      Flist[[i]]], 
AppendTo[babyFlist, 
                      Flist[[i]]]
], 
{b, 1, 
                    Nbbabies}
]
]; 
If[
                Length[ymates] ==  1  &&  
                  Length[bmates] ==  
                    0, 
Nybabies  =  
                  Random[PoissonDistribution[
                      aveybabies]]; 
Do[
If[
                    Random[] < 1/2, 
AppendTo[babyMylist, 
                      Flist[[i]]], 
AppendTo[babyFlist, 
                      Flist[[i]]]
], 
{b, 1, 
                    Nybabies}
]
]
], 
{i, 1, 
              Length[Flist]}
]; 
{babyMolist, babyMblist, babyMylist, 
            babyFlist}
];

Evolving lizards

\(Mo0\  = \ 30; \  (*\(70;\)*) 

    Molist\  = 
Table[
{Random[], 
          Random[]}, 
{i, 1, 
          Mo0}
];\)
Mb0  = 
      120;   (*140*) 
Mblist  = 
Table[
{Random[], Random[]}, 
{i, 1, 
          Mb0}
];
My0  =  
      30;   (*95;*) 
Mylist  = 
Table[
{Random[], 
          Random[]}, 
{i, 1, 
          My0}
];
F0  =  Mo0  +  Mb0  +  My0;
Flist  = 
Table[
{Random[], 
            Random[]}, 
{i, 1, 
            F0}
];
babyMolist  = {};
babyMblist  = {};
babyMylist  = {};
babyFlist  = {};


Show[
draw[
          Molist, {1, 0.25, 0}], 
draw[
          Mblist, {0, 0, 1}], 
draw[
          Mylist, {1, 1, 0}], 
draw[
          Flist, {0, 0, 0}], 
Prolog -> 
          PointSize[
            0.015], 
PlotRange -> {{-0.01, 
              1.01}, {-0.01, 1.01}}, 
Axes -> 
          False, 
AspectRatio -> 
          Automatic, 
DisplayFunction -> $DisplayFunction
];


Y  =  10;
y  =  0;
counts  =  {{0, Mo0, Mb0, My0, F0}};
While[
        y ≤ 
            Y  &&  (y == 0  ||  
              Last[Last[counts]] > 
                0), 
Do[
{Molist, 
              Mblist, Mylist, Flist, babyMolist, babyMblist, babyMylist, 
              babyFlist} =  
            move[Molist, Mblist, Mylist, Flist, babyMolist, babyMblist, 
              babyMylist, babyFlist], 
{w, 1, 
            26}
]; 
Do[
{newbabyMolist, newbabyMblist, 
              newbabyMylist, newbabyFlist}  = 
mate[
              Molist, Mblist, Mylist, 
              Flist]; 
babyMolist  =  
            Join[babyMolist, 
              newbabyMolist]; 
babyMblist  =  
            Join[babyMblist, 
              newbabyMblist]; 
babyMylist  =  
            Join[babyMylist, 
              newbabyMylist]; 
babyFlist  =  
            Join[babyFlist, newbabyFlist]; 
{Molist, 
              Mblist, Mylist, Flist, babyMolist, babyMblist, babyMylist, 
              babyFlist} =  
            move[Molist, Mblist, Mylist, Flist, babyMolist, babyMblist, 
              babyMylist, babyFlist], 
{w, 1, 
            26}
]; 
y  =  
          y + 1; 
Molist  =  
          Join[Molist, babyMolist]; 
Mblist  =  
          Join[Mblist, babyMblist]; 
Mylist  =  
          Join[Mylist, babyMylist]; 
Flist  =  
          Join[Flist, 
            babyFlist]; 
AppendTo[
counts, 
{y, 
            Length[Molist], Length[Mblist], Length[Mylist], 
            Length[Flist]}
]; 
babyMolist  =  {}; 
babyMblist  =  {}; 
babyMylist  =  {}; 
babyFlist  =  {}
];

Plotting the population of lizards

point  =  
        Graphics[{PointSize[0.00001], 
            Point[{0, 0}]}];
triangle  = 
Graphics[
{
Line[{{0, 0}, {2/Sqrt[3], 0}, {1/Sqrt[3], 1}, {0, 
                  0}}], 
{RGBColor[1, 0.25, 
                0], 
Point[{0, 
                  0}]
}, 
{RGBColor[0, 0, 
                1], 
Point[{1/Sqrt[3], 
                  1}]
}, 
{RGBColor[1, 1, 
                0], 
Point[{2/Sqrt[3], 
                  0}]
}
}
];


orbit  =  
Table[
{
(counts[[i, 3]] + 
                  2*counts[[i, 4]])/
(Sqrt[
                    3]*(counts[[i, 2]] + counts[[i, 3]] + 
                      counts[[i, 
                          4]])) , 
counts[[i, 
                  3]]/
(counts[[i, 2]] + 
                  counts[[i, 3]] + 
                  counts[[i, 
                      4]])
}, 
{i, 1, 
            Length[counts]}
];

orbitplot  = 
ListPlot[
orbit, 
Prolog -> {PointSize[0.015], 
              Thickness[0.005]}, 
PlotJoined -> 
            True, 
PlotStyle -> 
            RGBColor[1, 0, 0], 
DisplayFunction -> 
            Identity
];
Show[
triangle, 
        orbitplot, 
Prolog -> {PointSize[0.015], 
            Thickness[0.005]}, 
Axes -> 
          False, 
AspectRatio -> 
          Automatic, 
PlotRange -> 
          All, 
DisplayFunction -> $DisplayFunction
];

counts