ブナの成り年
ブナはドングリをつける樹木です。東北地方のブナは年によってたくさん実をつけたり、全く実をつけなかったりします。すなわち、翌年以降の数年間の結実が資源的に厳しくなるほどにある年の結実に非常に大きい投資をし、集団で一斉に開花します。樹木のなかではブナのように集団で揃って開花結実するものは比較的たくさんありますが、このようにするのが適応的な理由は受粉の効率が一斉開花によって高まるためであるとする説があります。今回はこのモデルに従ってシミュレーションをしていきます。このページの最後にC++のソースコードを掲載していますので、参考にしてもらえればと思います。
モデルの概要
まず、森林の中にN個の樹木が育つことのできる区画を想定します。そこに1本ずつ樹木が育つのですが、樹木は栄養状態の指標Yをそれぞれ持っています。Yが正の値であるとき結実できるとします。毎年、光合成による同化で1年あたり1ずつこの値が蓄積されます。幼木のうちは実をつけない設定ですが、成木はYが正であれば実をつけ、結実分だけYが減少します。結実量はYの値と他家花粉量に依存するようになっています。このようにしてモデル化された集団を交配と突然変異によって進化させていきます。進化はYの余剰分1に対してどれだけの種子を生産するかというパラメータkの変化によって観察できます。kが大きいと結実年と非結実年が顕著に現れるようになります。反対に、kが１に近いと毎年少しづつ結実するという性質を持つことを意味します。一斉開花結実が進化するかしないかは、パラメーターとして設定される、受粉成功の他家花粉量依存強度Beta、種子(実生)の生存率Sなどによって決定されると言われています。今回のシミュレーションモデルは、数理生物に関するYuuya Tachiki、Yoh Iwasaによる論文Both seedling banks and specialist seed predators promote the evolution of synchronized and intermittent reproduction (masting) in trees (Abstract)に沿って作成しています。詳しくはこの論文を参考にしてください。
受粉成功の他家花粉量依存強度が大きく、種子の生存率が高いと一斉開花結実は進化する
Beta=2.0、S=0.8というようなパラメーターを設定すると各個体が揃って開花結実するようになります。このような条件だとkが世代を経るごとに増加していきます。
受粉成功の他家花粉量依存強度が小さいと一斉開花結実は進化しない
Betaが小さいとき、一斉開花結実は進化しません。なぜなら、他個体の花粉が必要なければ自分だけで種子を作ることができるため、他個体と同年に集中して開花する必要がないからです。kは先の例よりも小さい値に収束しました。
種子の寿命を短くすると一斉開花結実は進化しない
Sがほとんど0に近いとき、一斉開花結実は進化しません。なぜなら、他個体と同調することで花粉の多さによって多くの種子を作ることができても、その年にできた種子の数は他個体のものも同様に豊作となるため、競争が激しくなる分だけメリットを打ち消してしまうからです。
ソースコード
このソースコードはGNU GPLライセンスにより保護されています。GNU GPLライセンスの条項に関してはこちらをご覧ください。
乱数にメルセンヌツイスタを採用しています。メルセンヌツイスタのソースコードはこちらを使用しています。このソースコード内ではヘッダーファイルMT.hとしてインクルードしています。ヘッダーファイルとして配布している場所もあります。
//
// masting_simulation
//
// Copyright 2015 Yuki Mochizuki
//
#define Year 10000000
#define N 100
#define Tau 10
#define Delta 0.04
#define S 0.0
#define Beta 2.0
#define M 0.01
#define M_width 0.1
#define MeanSeedPerYear 1
#include <iostream>
#include <math.h>
#include <vector>
#include "MT.h"
using namespace std;
class seed
{
public:
double k;
int age;
seed(){
k=1;
age=0;
}
};
int main() {
char *output ="masting_simulation.txt";
FILE* fpo;
fpo=std::fopen( output, "w" );
fprintf(fpo, "N %d, tau %d, Delta %f, S %f, M %f, Mw %f, Beta %f",N,Tau,Delta,S,M,M_width,Beta);
init_genrand(789749);
//data storage
int theta[N];
double k[N];
double Y[N];
int Yp[N];
double P[N];
int age[N];
vector<seed> seeds;
vector<seed> tmp_seeds;
vector<double> pollen_pool;
//temporal variables
int i=0,j=0,pollen=0,recuruit=0;
int seednum[N];
int seed_num=0,pollen_num=0,Yp_count=0;
double k_mean=0,age_mean=0;
seed seed;
//initiate
for(i=0;i<N;i++){
theta[i]=2;
Y[i]=(genrand_real2()-0.5)*3;
k[i]=genrand_real2()+0.5;
age[i]=(int)(genrand_real2()*100.0);
if(age[i]>Tau)theta[i]=2;
else theta[i]=1;
}
for(int y=0;y<Year;y++){
//growth phase
for(i=0;i<N;i++){
Y[i]++;
}
//reproduction phase
pollen=0;
for(i=0;i<N;i++){
if(Y[i]>0)Yp[i]=Y[i];
else Yp[i]=0;
}
for(i=0;i<N;i++){
pollen+=(theta[i]==2)*Yp[i];
}
for(i=0;i<N;i++){
P[i]=pow((double)(pollen-Yp[i])/(double)(N-1),Beta);
}
vector<double>().swap(pollen_pool);
for(i=0;i<N;i++){
seednum[i]=k[i]*P[i]*Yp[i];
for(j=0;j<seednum[i];j++){
pollen_pool.push_back(k[i]);
}
}
pollen_num=pollen_pool.size();
for(i=0;i<N;i++){
Y[i]-=seednum[i];
for(j=0;j<seednum[i]*MeanSeedPerYear;j++){
if(genrand_real2()<M){
seed.k=0.5*(k[i]+pollen_pool.at((int)(genrand_real2()*pollen_num)))+(genrand_real2()-0.5)*2*M_width;
seeds.push_back(seed);
}else{
seed.k=0.5*(k[i]+pollen_pool.at((int)(genrand_real2()*pollen_num)));
seeds.push_back(seed);
}
}
}
seed_num=(int)seeds.size();
//death and recruitment phase
for(i=0;i<N;i++){
if(theta[i]==0||genrand_real2()<Delta){
if(seed_num>0){
recuruit=(int)((double)seed_num*genrand_real2());
k[i]=seeds.at(recuruit).k;
Y[i]=0;
age[i]=0;
theta[i]=1;
}else{
theta[i]=0;
}
}
}
//seed death phase
if(seed_num>0){
vector<class seed>().swap(tmp_seeds);
for(i=0;i<seed_num;i++){
if(genrand_real2()<S){
tmp_seeds.push_back(seeds.at(i));
}
}
vector<class seed>().swap(seeds);
seeds=tmp_seeds;
}
//age
for(i=0;i<N;i++){
if(theta[i]>0)age[i]++;
if(age[i]>Tau)theta[i]=2;
}
//output
if(y%1000==0){
k_mean=0;
age_mean=0;
Yp_count=0;
for(i=0;i<N;i++){
if(theta[i]>0){
k_mean+=k[i];
age_mean+=age[i];
if(Y[i]>0)Yp_count++;
}
}
k_mean/=(double)N;
age_mean/=(double)N;
cout << "Year=" << y << " k(mean)=" << k_mean << " age(mean)=" << age_mean << " seeds=" << seed_num << " fruit=" << Yp_count<< endl;
fprintf(fpo, "%d %f %f %d %d\n",y,k_mean,age_mean,seed_num,Yp_count);
fflush(fpo);
}
}
return 0;
}