Immigration - emigration model

We first introduce a deterministic dynamics given as ODE and then extend it to a stochastic dynamics.

Immigration model with a constant inflow

We assume a population to which individuals immigrate (flow in) with a constant rate.

Deterministic model as ODE: pop size n(t) is a real

deq={n'[t]==flowIn,n[0]==n0}

sol=DSolve[deq,n[t],t]

para={flowIn->10.0,n0->10}

deqN=deq/.para

solN=NDSolve[deqN,n[t],{t,0,10}]

Plot[Evaluate[n[t]/.solN],{t,0,10},AxesLabel->{"t","n(t)"}]

This ODE is too easy and Mathematica would not be necessary. The solution is that the population size n(t) linearly increases with time t starting from an initial value n0.

Stochastic model: pop size n(t) is an integer

We now extend the deterministic ODE model as a stochastic model where population size is integer and population size (state variable) jumps from n to n + 1 stochastically.

The rule is that time to the next event (an individual joins the population) follows the exponential distirbution with the constant parameter.

rateInflow=2;

Table[-Log[Random[]]/rateInflow,{10}]

This sequence gives a series of inter-event times (waiting times to the next event). After each inter-event time, population size is incremented by one. The stochastic process can be implemented as follows.

rateInflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)

t=0;  (* Initial condition *)
popSize=0;
seq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateInflow;
t+=wt;
popSize++;
AppendTo[seq,{t,popSize}]
]

seq

ListPlot[seq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Note that the population size increases by one at each time an inflow occurs. Because this is a stochastic process, each realization results in different dynamics. 
To view the stochastic dynamics, it is often convenient to observe the population size with a FIXED time interval, say DT = 1.0. To do this, we introduce a time tEq that always advances by a constant DT and record the population size at each tEq.

rateInflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=0;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateInflow;

(* While loop to advance tEq with DT *)
While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
popSize++;
AppendTo[seq,{t,popSize}]
]

seq

seqEq

Let' s compare the two sequence data, seq and seqEq. Both should be the same but they differ in the timing of population census; “seq” is original record of population change. In “seqEq”, we census population size with a FIXED time interval DT, i.e., t = 0, DT, 2DT, 3Dt, ...

ListPlot[{seq,seqEq},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

If we run the process many times with the same conditions, we may have a hint to understand the stochastic process.Let' s run the simulation multiple times.

rateInflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* Number of realizations *)

seqSims={};  (* Sequence that records the population size in each realization *)
seqSimsEq={}; (* Sequence that records the population size measured at tEq in each realization *)

Do[
t=0;
tEq=0;
popSize=0;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateInflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
popSize++;
AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSims

seqSimsEq

ListPlot[seqSims,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Take out the time sequence

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]]

Take out the sequences of the population size

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq]

Transpose[seqPopSize]

Calculate the average (or mean) of the population size censusced at tEq = 0, 1, 2, 3, ...

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]]

Combine the time sequence with the sequence of the averaged population size

seqTPopSizeAverage=Transpose[{seqT,seqPopSizeAverage}]

gSimAverage=ListPlot[seqTPopSizeAverage]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeAverage,{1,t},t]

gFit=Plot[fit,{t,0,10}]

Show[gSimAverage,gFit]

How does the averaged population size change with time? 
Increase the number of realizations and check the result. 
Slope of the linear regression should converge to the “inflow rate” for a large number of realizations. For analytical proof, see the documents at http://gi.ics.nara-wu.ac.jp/~takasu/lecture/graduate20.html

Emigration model with a constant outflow

We assume a population from whicn individuals emigrate (flow out) with a constant rate.

Deterministic model as ODE: pop size n(t) is a real

deq={n'[t]==-flowOut,n[0]==n0}

sol=DSolve[deq,n[t],t]

para={flowOut->10.0,n0->100}

deqN=deq/.para

solN=NDSolve[deqN,n[t],{t,0,10}]

Plot[Evaluate[n[t]/.solN],{t,0,10},AxesLabel->{"t","n(t)"}]

This ODE is too easy and Mathematica would not be necessary. The solution is that the population size n(t) linearly increases with time t starting from an initial value n0.

Stochastic model: pop size n(t) is an integer

We now extend the deterministic ODE model as a stochastic model where population size is integer and population size (state variable) jumps from n to n - 1 stochastically.

The rule is that time to the next event (an individual joins the population) follows the exponential distirbution with the constant parameter.

rateOutflow=2;

Table[-Log[Random[]]/rateOutflow,{10}]

This sequence gives a series of inter-event times (waiting times to the next event). After each inter-event time, population size is incremented by one. The stochastic process can be implemented as follows.

rateOutflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)

t=0;  (* Initial condition *)
popSize=100;
seq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;
t+=wt;
popSize--;
AppendTo[seq,{t,popSize}]
]

seq

ListPlot[seq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Note that the population size increases by one at each time an inflow occurs. Because this is a stochastic process, each realization results in different dynamics. 
To view the stochastic dynamics, it is often convenient to observe the population size at a fixed time interval, say DT = 1.0. To do this, we introduce a time tEq that always advances by a constant DT and record the population size at each tEq.

rateOutflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=100;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
popSize--;
AppendTo[seq,{t,popSize}]
]

seq

seqEq

ListPlot[{seq,seqEq},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

If we run the process many times with the same conditions, we may have a hint to understand the stochastic process.Let' s run the simulation multiple times.

rateOutflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* Number of realizations *)

seqSims={};  (* Sequence that records the population size in each realization *)
seqSimsEq={}; (* Sequence that records the population size measured at tEq in each realization *)

Do[
t=0;
tEq=0;
popSize=100;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
popSize--;
AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSims

seqSimsEq

ListPlot[seqSims,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Take out the time sequence

seqSimsEq[[1]]

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]]

Take out the sequences of the population size

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq]

Transpose[seqPopSize]

Calculate the average (or mean) of the population size measured at tEq = 0, 1, 2, 3, ...

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]]

Combine the time sequence with the sequence of the averaged pouplation size

seqTPopSizeAverage=Transpose[{seqT,seqPopSizeAverage}]

gSimAverage=ListPlot[seqTPopSizeAverage]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeAverage,{1,t},t]

gFit=Plot[fit,{t,0,10}]

Show[gSimAverage,gFit]

How does the averaged population size change with time? 
Increase the number of realizations and check the result. 
Slope of the linear regression should converge to the “outflow rate” for a large number of realizations. For analytical proof, see the documents at http://gi.ics.nara-wu.ac.jp/~takasu/lecture/graduate20.html

Stochastic model: pop size n(t) is an integer but it cannot be negative

We now extend the deterministic ODE model as a stochastic model where population size is integer and population size (state variable) jumps from n to n - 1 stochastically.

The rule is that time to the next event (an individual joins the population) follows the exponential distribution with the constant parameter.

rateOutflow=2;

Table[-Log[Random[]]/rateOutflow,{10}]

This sequence gives a series of inter-event times (waiting times to the next event). After each inter-event time, population size is incremented by one. The stochastic process can be implemented as follows.

rateOutflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)

t=0;  (* Initial condition *)
popSize=10;
seq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;
t+=wt;
popSize--;

If[popSize < 0,popSize = 0];

AppendTo[seq,{t,popSize}]
]

seq

ListPlot[seq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Note that the population size increases by one at each time an inflow occurs. Because this is a stochastic process, each realization results in different dynamics. 
To view the stochastic dynamics, it is often convenient to observe the population size at a fixed time interval, say DT = 1.0. To do this, we introduce a time tEq that always advances by a constant DT and record the population size at each tEq.

rateOutflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=10;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
popSize--;

If[popSize < 0,popSize = 0];

AppendTo[seq,{t,popSize}]
]

seq

seqEq

ListPlot[{seq,seqEq},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

If we run the process many times with the same conditions, we may have a hint to understand the stochastic process.Let' s run the simulation multiple times.

rateOutflow = 2;  (* Set the inflow rate *)
tEnd=10;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* Number of realizations *)

seqSims={};  (* Sequence that records the population size in each realization *)
seqSimsEq={}; (* Sequence that records the population size measured at tEq in each realization *)

Do[
t=0;
tEq=0;
popSize=10;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
popSize--;

If[popSize < 0,popSize = 0];

AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSims

seqSimsEq

ListPlot[seqSims,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Take out the time sequence

seqSimsEq[[1]]

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]]

Take out the sequences of the population size

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq]

Transpose[seqPopSize]

Calculate the average (or mean) of the population size measured at tEq = 0, 1, 2, 3, ...

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]]

Combine the time sequence with the sequence of the averaged population size

seqTPopSizeAverage=Transpose[{seqT,seqPopSizeAverage}]

gSimAverage=ListPlot[seqTPopSizeAverage]

Immigration and emigration model

We assume a population with immigration and emigration.

Deterministic model as ODE: pop size n(t) is a real

deq={n'[t]==flowIn-flowOut,n[0]==n0}

sol=DSolve[deq,n[t],t]

para={flowIn->12,flowOut->10,n0->100}

deqN=deq/.para

solN=NDSolve[deqN,n[t],{t,0,10}]

Plot[Evaluate[n[t]/.solN],{t,0,10},AxesLabel->{"t","n(t)"}]

This ODE is too easy and Mathematica would not be necessary. The solution is that the population size n(t) linearly increases with time t starting from an initial value n0.

Stochastic model: pop size n(t) is an integer: Simulation and the moment dynamics

We now extend the deterministic ODE model as a stochastic model where population size is integer and population size (state variable) jumps from n to n + 1 with constant inflow and jumps from n to n – 1 with constant outflow stochastically.

The rule is that time to the next event, an individual is coming in or out, follows the exponential distribution with the constant “inflow + outflow”.

And after the inter-event time, we have to choose which event, an individual coming in or out, has occurred. Here we use the rule that immigration (in) and emigration (out) occurs with the probability proportional to the inflow and the outflow, respectively.

rateInflow = 2.5;
rateOutflow=2;

Table[-Log[Random[]]/(rateInflow+rateOutflow),{10}]

This sequence gives a series of inter-event times (waiting times to the next event). After each inter-event time, population size is incremented by one with the probability rateInflow/(rateInflow + rateOutflow) and is decreased by one with the rest of the probability. The stochastic process can be implemented as follows.

rateInflow = 3;
rateOutflow = 2;  (* Set the inflow rate *)
tEnd=50;        (* Set the end of time *)

t=0;  (* Initial condition *)
popSize=10;
seq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateInflow+rateOutflow);
t+=wt;
If[Random[]<rateInflow/(rateInflow+rateOutflow),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];
AppendTo[seq,{t,popSize}]
]

seq;

ListPlot[seq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Note that the population size increases by one at each time an inflow occurs. Because this is a stochastic process, each realization results in different dynamics. 
To view the stochastic dynamics, it is often convenient to observe the population size with a FIXED time interval, say DT = 1.0. To do this, we introduce a time tEq that always advances by a constant DT and record the population size at each tEq.

rateInflow = 3; (* Set the inflow rate *)
rateOutflow = 2;  (* Set the outflow rate *)
tEnd=50;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=10;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateInflow+rateOutflow);

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;
If[Random[]<rateInflow/(rateInflow+rateOutflow),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];
AppendTo[seq,{t,popSize}]
]

seq;

seqEq

ListPlot[{seq,seqEq},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

If we run the process many times with the same conditions, we may have a hint to understand the stochastic process.Let' s run the simulation multiple times.

rateInflow = 3; (* Set the inflow rate *)
rateOutflow = 2;  (* Set the outflow rate *)
tEnd=50;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=100;  (* Number of realizations *)

seqSims={};  (* Sequence that records the population size in each realization *)
seqSimsEq={}; (* Sequence that records the population size measured at tEq in each realization *)

Do[
t=0;
tEq=0;
popSize=10;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateInflow+rateOutflow);

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;

If[Random[]<rateInflow/(rateInflow+rateOutflow),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];

AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSims

seqSimsEq;

ListPlot[seqSims,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Take out the time sequence

seqSimsEq[[1]]

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]]

Take out the sequences of the population size

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq];

Transpose[seqPopSize];

Calculate the average (or mean) of the population size measured at tEq = 0, 1, 2, 3, ...

seqPopSizeAverage=Map[Mean,Transpose[seqPopSize]]

Combine the time sequence with the sequence of the averaged population size

seqTPopSizeAverage=Transpose[{seqT,seqPopSizeAverage}]

gSimAverage=ListPlot[seqTPopSizeAverage]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeAverage,{1,t},t]

gFit=Plot[fit,{t,0,50}]

Show[gSimAverage,gFit]

seqPopSizeVar=Map[Variance,Transpose[seqPopSize]]

Combine the time sequence with the sequence of the averaged population size

seqTPopSizeVar=Transpose[{seqT,seqPopSizeVar}]

gSimVar=ListPlot[seqTPopSizeVar]

Erase the variable t

t=.

Linear regression of the averaged population size

fit=Fit[seqTPopSizeVar,{1,t},t]

gFit=Plot[fit,{t,0,50}]

Show[gSimVar,gFit]

How does the averaged population size change with time? 
Increase the number of realizations and check the result. 
Slope of the linear regression should converge to the “inflow rate” - “outflow rate” for a large number of realizations. For analytical proof, see the documents at http://gi.ics.nara-wu.ac.jp/~takasu/lecture/graduate20.html

Stochastic model: pop size n(t) is an integer but it cannot be negative

We now extend the deterministic ODE model as a stochastic model where population size is integer and population size (state variable) jumps from n to n + 1 with constant inflow and jumps from n to n – 1 with constant outflow stochastically.

The rule is that time to the next event, an individual is coming in or out, follows the exponential distribution with the constant “inflow + outflow”.

And after the inter-event time, we have to choose which event, an individual coming in or out, has occurred. Here we use the rule that immigration (in) and emigration (out) occurs with the probability proportional to the inflow and the outflow, respectively.

We also restict the population size to be non-negative. This makes a biological sense because no individual can get out of empty population.

This sequence gives a series of inter-event times (waiting times to the next event). After each inter-event time, population size is incremented by one with the probability rateInflow/(rateInflow + rateOutflow) and is decreased by one with the rest of the probability. The stochastic process can be implemented as follows.

Here we restrict the population size to be non-negative.

rateInflow = 2;
rateOutflow = 3;  (* Set the inflow rate *)
tEnd=20;        (* Set the end of time *)

t=0;  (* Initial condition *)
popSize=10;
seq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/(rateInflow+rateOutflow);
t+=wt;
If[Random[]<rateInflow/(rateInflow+rateOutflow),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];

If[popSize<0,popSize=0];  (* We restrict popSize non-nagative *)

AppendTo[seq,{t,popSize}]
]

seq

ListPlot[seq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Note that the population size increases by one at each time an inflow occurs. Because this is a stochastic process, each realization results in different dynamics. 
To view the stochastic dynamics, it is often convenient to observe the population size at a fixed time interval, say DT = 1.0. To do this, we introduce a time tEq that always advances by a constant DT and record the population size at each tEq.

rateInflow = 2; (* Set the inflow rate *)
rateOutflow = 3;  (* Set the outflow rate *)
tEnd=20;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

t=0;
tEq=0;  (* Time that advances by a constant DT *)
popSize=10;
seq={};
seqEq={}; (* Sequence that records the population size measured at tEq *)

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;

If[Random[]<rateInflow/(rateInflow+rateOutflow),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];

If[popSize<0,popSize=0];(* We restrict popSize non-nagative *)

AppendTo[seq,{t,popSize}]
]

seq

seqEq

ListPlot[{seq,seqEq},Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

If we run the process many times with the same conditions, we may have a hint to understand the stochastic process.Let' s run the simulation multiple times.

rateInflow = 2; (* Set the inflow rate *)
rateOutflow = 3;  (* Set the outflow rate *)
tEnd=30;        (* Set the end of time *)
DT=1;     (* Constant time interval *)

numRealizations=200;  (* Number of realizations *)

seqSims={};  (* Sequence that records the population size in each realization *)
seqSimsEq={}; (* Sequence that records the population size measured at tEq in each realization *)

Do[
t=0;
tEq=0;
popSize=10;
seq={};
seqEq={};

AppendTo[seq,{t,popSize}];

While[t<tEnd,
wt=-Log[Random[]]/rateOutflow;

While[ t+wt > tEq,
AppendTo[seqEq,{tEq,popSize}];
tEq+=DT;
If[tEq>tEnd,Break[]];
];

t+=wt;

If[Random[]<rateInflow/(rateInflow+rateOutflow),popSize++,  (* Here we choose inflow or outflow stochastically *)
popSize--
];

If[popSize<0,popSize=0];(* We restrict popSize non-nagative *)

AppendTo[seq,{t,popSize}]
];
AppendTo[seqSims,seq];
AppendTo[seqSimsEq,seqEq],{i,numRealizations}
]

seqSims

seqSimsEq

ListPlot[seqSims,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

ListPlot[seqSimsEq,Joined->True,PlotMarkers->Automatic,AxesLabel->{"t","PopSize"}]

Take out the time sequence

seqSimsEq[[1]]

seqT=Map[Transpose[#][[1]]&,seqSimsEq][[1]]

Take out the sequences of the population size

seqPopSize=Map[Transpose[#][[2]]&,seqSimsEq]

popSizeLast=Last[Transpose[seqPopSize]]

Max[popSizeLast]

gHist=Histogram[popSizeLast,{0,15,1},"PDF"]

gDist=DiscretePlot[PDF[GeometricDistribution[1/3],k],{k,0,15},PlotMarkers->Automatic,PlotRange->All]

Show[gHist,gDist]

When rateInflow < rateOutflow and the population size is restricted to be non-negative, population size after a long run shows a certain distribution.

It can be shown the distribution is the geometric (exponential) distribution with the parameter (rateOutflow – rateInflow)/rateOutflow. 

For analytical proof, see the documents at http://gi.ics.nara-wu.ac.jp/~takasu/lecture/graduate20.html

This immigration-emigration model is essentially the same as a model in “Queueing theory”. A familiar example is a cash register at supermarket where customers arrive at the register with a constant rate and the register handles each customer with a constant handling rate. It’s obvious that the handling rate should be larger than the arrival rate. Otherwise, waiting cue at the register grows and customers will complain.

QueueingProcess function of Mathematica

?QueueingProcess

data=RandomFunction[QueueingProcess[3,3],{0,100}]

ListStepPlot[data,Filling->Axis]

data["LastValues"]

lastValuesSeq={};
Do[
data=RandomFunction[QueueingProcess[2,3],{0,100}];
AppendTo[lastValuesSeq,data["LastValues"][[1]]],{200}
]

lastValuesSeq

Max[lastValuesSeq]

Histogram[lastValuesSeq,{0,20,1},"PDF"]