myappinstaller_web.exe |
When YY Brook Trout Impact Native Bull Trout
Matt Falcy
November 8, 2017
Background
Non-native brook trout can displace native bull trout. Efforts to control brook trout can inlcude the release of YY male brook trout. Such individuals will produce all male offspring, thereby skewing the sex ratio of the brook trout population and supressing population size. However, since brook trout and bull trout are in close proximity and readily hybridize, there is concern that introducing YY brook trout can adversely affect bull trout.
My objective to is to identify "conditions" where release of YY brook trout can be problematic for bull trout. Two kinds of models are developed below: a system of coupled ordinary differential equations (CODE) and an individual-based model (IBM). I will show that these models can produce idential results as a "proof" of the IBM, then extend the IBM in a way that would be prohibitive using CODE.
The System
Let F and M denote females and males, respectively. Let x and y denote the autosomal karyotope of individuals. For natural salmoninds, Fxx and Mxy mating will produce 50% Fxx and 50% Mxy. However, when artificially produced Myy breed with Fxx then all offspring are Mxy.
Both the CODE and IBM will track these three kinds of individuals (sexs). Since individuals readily hybridize, I will not track individual geneologies. Rather, there will be three locations. These locations should be thought of as a brook trout zone, a bull trout zone, and a hybrid zone. However, movement among these zones and subsequent breeding will result in all three zones containing genetic legacy from both bull and brook trout.
The important response variable to keep track of is the abudnance of females in the putative "bull" zone.
CODE- Coupled Ordinary Differential Equations
Females are always karyotpe xx, so simply denote them as "F". Use subscripts 1, 2, 3 to denote zones. I is a 3 by 3 matrix giving movement probabilities between the three zones. Note that all individuals in a zone, regardless of sex, have the same movement probablities (for now). is a birth coefficient in zone i, and is a death coefficient in zone i. is the carrying capacity in zone i. %\mu_i% is the number of yy males introduced into zone i by humans.
The modeling presented below is similar to Gutierrez and Teem (2006) Journal of Thoeretical Biology, 241:333-341.
<include>metapop.m</include>
clear t y
K=[300 500 100];
D=[0.08 0.1 0.12];
B=[0.08 0.05 0.05];
mu=5;
I=[NaN 0 0; %Disallow movement (for now)
0 NaN 0;
0 0 NaN];
Tmax=100;
[t,y]=ode45(@(t,y) metapop(t,y,K,B,D,mu,I),[0 Tmax],[25,25,10,25,25,10,10,10,10]);
subplot(1,2,1),plot(t,y(:,2))
hold on
subplot(1,2,1),plot(t,y(:,5))
subplot(1,2,1),plot(t,y(:,8))
xlim([0 Tmax])
title({'Coupled Differential Equations', 'no movement'})
xlabel('Time t');
ylabel('Female Abundance');
legend('Brook','Bull','Hybrid')
IBM- Individual-Based Model
% Initialize metapopulation
clear all
N0=[25 25 10;
25 25 10;
10 10 10];%rows are sex, cols are location
N=[];
for i=1:3 %sex
for j=1:3 %location
N((size(N,1)+1):(size(N,1)+N0(i,j)),1:2)=repmat([i,j],N0(i,j),1);
end
end
%{
%stepping stone
M(:,:,1)=[0.70, 0, 0.3;
0, 1, 0.0;
0, 0.3, 0.7];%females
M(:,:,2)=[0.70, 0, 0.3;
0, 1, 0.0;
0, 0.3, 0.7];%males
M(:,:,3)=[0.70, 0, 0.3;
0, 1, 0.0;
0, 0.3, 0.7];%supermales
%}
%Disallow movement as a null scenario
M(:,:,1)=diag([1,1,1]);
M(:,:,2)=diag([1,1,1]);
M(:,:,3)=diag([1,1,1]);
%Make sure rows sum to 1
for i=1:3
sum(M(:,:,i),2)
end
Tmax=100;
out1=repmat(NaN,[3,3,Tmax]);
out2=repmat(NaN,[3,3,Tmax]);
for t=1:Tmax
% Movement
for i=1:length(N)
pd=makedist('Multinomial','Probabilities',M(N(i,2),:,N(i,1)));
N(i,2)=random(pd);
end
%Death
% Set death rate as a constant and apply density effects on reproduction as per
% Gutierrez and Teem (2006) Journal of Thoeretical Biology, 241:333-341.
% Later, we should place density effects on deat rate. Later, we can specifiy unique
% death rates depending on sex.
% Set death rate in each area(Brook, Bull, Hybrid)
D=[0.08, 0.1, 0.12];
clear dead
for i=1:length(N)
dead(i)=binornd(1,D(N(i,2)));
end
N=N(~dead,:);%if not dead then alive.
out1(:,:,t)=crosstab(N(:,1),N(:,2));%rows are sexes, cols are areas.
%Reproduction
%Sum all individuals in each area
%Ni=hist(N(:,2),3);
Ni=sum(crosstab(N(:,1),N(:,2)));
if (sum(Ni)==length(N))==0, disp('Totals dont add up')
end
% Set carrying capacity in each area (Brook, Bull, Hybrid)
K=[300, 500, 100];
% Compute simple logistic term. G&T(2006)
L=1-Ni./K;
%replicate G&T
B=[0.08 0.05 0.05];% Set birth rate in each area (Brook, Bull, Hybrid).
Z=crosstab(N(:,1),N(:,2));%rows are sexes, cols are areas.
Supers=[5,0,0]; %only place supermales in Zone 1
for i =1:3 %cycle over locations
N_f(i)=floor(0.5*B(i)*Z(1,i)*Z(2,i)*L(i));%Females in area i. Half the matings (0.5) are female.
N_m(i)=floor((0.5*B(i)*Z(1,i)*Z(2,i)+B(i)*Z(1,i)*Z(3,i))*L(i));%Males in area i. Note that supermales produce all males.
N_s(i)=Supers(i); %Constant dump of supermales by humans.
end
%If below capacity
for i=1:3 %cycle over locations
N=vertcat(N,repmat([1,i],[N_f(i),1]));%Females
N=vertcat(N,repmat([2,i],[N_m(i),1]));%Males
N=vertcat(N,repmat([3,i],[N_s(i),1]));%Supermales
end
%{
%If above capacity, implement 'negative reproduction'.
for i=1:3 %cycle over areas
if N_f(i)<0 %females
idx=find(N(:,1)==1 & N(:,2)==i);
N(idx(1:abs(N_f(i))),:)=[];%Negative reproduction
end
if N_m(i)<0 %males
idx=find(N(:,1)==2 & N(:,2)==i);
N(idx(1:abs(N_m(i))),:)=[];%Negative reproduction
end
if N_s(i)<0 %supermales
idx=find(N(:,1)==3 & N(:,3)==i);
N(idx(1:abs(N_s(i))),:)=[];%Negative reproduction
end
end
%}
out2(:,:,t)=crosstab(N(:,1),N(:,2));%rows are sexes, cols are areas.
end %t loop
%{
figure
plot(1:Tmax,reshape(out1(1,1,:),Tmax,1))
hold on
plot(1:Tmax,reshape(out1(1,2,:),Tmax,1))
plot(1:Tmax,reshape(out1(1,3,:),Tmax,1))
xlabel('Time')
ylabel('Female Abundance')
legend('Zone 1','Zone 2','Zone 3')
%}
%
subplot(1,2,2),plot(1:Tmax,reshape(out1(1,1,:),Tmax,1))
hold on
subplot(1,2,2),plot(1:Tmax,reshape(out1(1,2,:),Tmax,1))
subplot(1,2,2),plot(1:Tmax,reshape(out1(1,3,:),Tmax,1))
xlim([0 Tmax])
title({'Individual Based Model', 'no movement'})
xlabel('Time')
ylabel('Female Abundance')
legend('Brook','Bull','Hybrid')
%}
Commentary
The IBM and CODE seem to prododuce qualitatively similar results. Only the IBM has stochasticity.
Note that addition of yy males into the brook zone suppresses the female abundance there to less than 50% of the carrying capacity (K=300). YY males are not introduced into the brook and hybrid zones, and so females abudnance approaches 50% of carrying capacity.
Now we will repeat the modeling presented above, but allowing fish to move (1) from brook zone to the hybrid zone, and (2) from the hybrid zone to the bull zone. Thus fish do NOT move (1) from brook to bull zones, (2) from hybrid to brook zones, (3) from the bull zone to either brook or hybrid zones. We will set the movement probability at 30%, and apply it uniformly to all sexes within the zone. Call this "Stepping Stone Movement" pattern.
Stepping Stone Movement
See the commentary immediately above for the movement rules.
clear t y
K=[300 500 100];
D=[0.08 0.1 0.12];
B=[0.08 0.05 0.05];
mu=5;
I=[NaN 0 0; %Stepping Stone
0 NaN 0.3;
0.3 0 NaN];
Tmax=100;
[t,y]=ode45(@(t,y) metapop(t,y,K,B,D,mu,I),[0 Tmax],[25,25,10,25,25,10,10,10,10]);
%
% <include>metapop.m</include>
figure
subplot(1,2,1),plot(t,y(:,2))
hold on
subplot(1,2,1),plot(t,y(:,5))
subplot(1,2,1),plot(t,y(:,8))
xlim([0 Tmax])
title({'Coupled Differential Equations', 'stepping stone movement'})
xlabel('Time t');
ylabel('Female Abundance');
legend('Brook','Bull','Hybrid')
% IBM- Individual-Based Model
% Initialize metapopulation
clear all
N0=[25 25 10;
25 25 10;
10 10 10];%rows are sex, cols are location
N=[];
for i=1:3 %sex
for j=1:3 %location
N((size(N,1)+1):(size(N,1)+N0(i,j)),1:2)=repmat([i,j],N0(i,j),1);
end
end
%stepping stone
M(:,:,1)=[0.70, 0, 0.3;
0, 1, 0.0;
0, 0.3, 0.7];%females
M(:,:,2)=[0.70, 0, 0.3;
0, 1, 0.0;
0, 0.3, 0.7];%males
M(:,:,3)=[0.70, 0, 0.3;
0, 1, 0.0;
0, 0.3, 0.7];%supermales
%{
%Disallow movement as a null scenario
M(:,:,1)=diag([1,1,1]);
M(:,:,2)=diag([1,1,1]);
M(:,:,3)=diag([1,1,1]);
%}
%Make sure rows sum to 1
for i=1:3
sum(M(:,:,i),2)
end
Tmax=100;
out1=repmat(NaN,[3,3,Tmax]);
out2=repmat(NaN,[3,3,Tmax]);
for t=1:Tmax
% Make them move probabilistically
for i=1:length(N)
pd=makedist('Multinomial','Probabilities',M(N(i,2),:,N(i,1)));
N(i,2)=random(pd);
end
%Death
% Set death rate as a constant and apply density effects on reproduction as per
% Gutierrez and Teem (2006) Journal of Thoeretical Biology, 241:333-341.
% Later, we should place density effects on deat rate. Later, we can specifiy unique
% death rates depending on sex.
% Set death rate in each area(Brook, Bull, Hybrid)
D=[0.08, 0.1, 0.12];
clear dead
for i=1:length(N)
dead(i)=binornd(1,D(N(i,2)));
end
N=N(~dead,:);%if not dead then alive.
out1(:,:,t)=crosstab(N(:,1),N(:,2));%rows are sexes, cols are areas.
%Reproduction
%Sum all individuals in each area
%Ni=hist(N(:,2),3);
Ni=sum(crosstab(N(:,1),N(:,2)));
if (sum(Ni)==length(N))==0, disp('Totals dont add up')
end
% Set carrying capacity in each area (Brook, Bull, Hybrid)
K=[300, 500, 100];
% Compute simple logistic term. G&T(2006)
L=1-Ni./K;
%replicate G&T
B=[0.08 0.05 0.05];% Set birth rate in each area (Brook, Bull, Hybrid).
Z=crosstab(N(:,1),N(:,2));%rows are sexes, cols are areas.
Supers=[5,0,0]; %only place supermales in Zone 1
for i =1:3 %cycle over locations
N_f(i)=floor(0.5*B(i)*Z(1,i)*Z(2,i)*L(i));%Females in area i. Half the matings (0.5) are female.
N_m(i)=floor((0.5*B(i)*Z(1,i)*Z(2,i)+B(i)*Z(1,i)*Z(3,i))*L(i));%Males in area i. Note that supermales produce all males.
N_s(i)=Supers(i); %Constant dump of supermales by humans.
end
%If below capacity
for i=1:3 %cycle over locations
N=vertcat(N,repmat([1,i],[N_f(i),1]));%Females
N=vertcat(N,repmat([2,i],[N_m(i),1]));%Males
N=vertcat(N,repmat([3,i],[N_s(i),1]));%Supermales
end
%{
%If above capacity, implement 'negative reproduction'.
for i=1:3 %cycle over areas
if N_f(i)<0 %females
idx=find(N(:,1)==1 & N(:,2)==i);
N(idx(1:abs(N_f(i))),:)=[];%Negative reproduction
end
if N_m(i)<0 %males
idx=find(N(:,1)==2 & N(:,2)==i);
N(idx(1:abs(N_m(i))),:)=[];%Negative reproduction
end
if N_s(i)<0 %supermales
idx=find(N(:,1)==3 & N(:,3)==i);
N(idx(1:abs(N_s(i))),:)=[];%Negative reproduction
end
end
%}
out2(:,:,t)=crosstab(N(:,1),N(:,2));%rows are sexes, cols are areas.
end %t loop
%{
figure
plot(1:Tmax,reshape(out1(1,1,:),Tmax,1))
hold on
plot(1:Tmax,reshape(out1(1,2,:),Tmax,1))
plot(1:Tmax,reshape(out1(1,3,:),Tmax,1))
xlabel('Time')
ylabel('Female Abundance')
legend('Zone 1','Zone 2','Zone 3')
%}
subplot(1,2,2),plot(1:Tmax,reshape(out2(1,1,:),Tmax,1))
hold on
subplot(1,2,2),plot(1:Tmax,reshape(out2(1,2,:),Tmax,1))
subplot(1,2,2),plot(1:Tmax,reshape(out2(1,3,:),Tmax,1))
xlim([0 Tmax])
title({'Individual Based Model', 'stepping stone movement'})
xlabel('Time')
ylabel('Female Abundance')
legend('Brook','Bull','Hybrid')