RでAdaptive Dynamics（工事中）

Adaptive Dynamics理論が扱う進化動態をR言語で楽にシミュレーションするサイトです。

• 様々な面白い進化的分岐を、短めのRプログラムで楽に速く計算・可視化することを重視します。
• ここで紹介するプログラムは、適応度関数の部分を変更すれば他のモデルの解析に使えるように書きます。自由に改変して自分のモデルの解析に使ってみて下さい。
• 大体のOSで楽にインストールできるのでRを選びました（Rのインストール解説（本家））。
• 多次元の色々な進化的分岐の計算プログラムを追加する予定です。

１次元形質空間の進化的分岐

• モデル：MacArthur-Levinsの資源競争モデル (例えばDieckmann and Doebeli 1999)
  
• プログラムファイル"branch1d3.R"(100行くらい)を保存して、そのディレクトリ（フォルダ）でRを起動して以下のコマンドを入力すれば計算が実行されグラフが表示されます（日本語のコメントが文字化けする場合は下に貼った中身をコピペして下さい）。windowsの場合はRを起動して「Rソースコードの読み込み」のような選択⇒branch1d3.Rでも実行できるはず。

       
  source("branch1d3.R")
  

• Adaptive Dynamics理論で解析すると、s_aがs_Kよりも小さい場合にx=0が進化的分岐点になり、計算結果を良く予測します（この解説の1.10節）。
  
• 表現型（個体）間の相互作用の計算が最も時間がかかるので以下のような工夫をしています。
  • 変異型の出現イベントは数百〜数千ステップおきにして、存在している表現型のみで個体数動態をまわす。
  • 相互作用計算においてexp関数を呼ぶ回数を可能な限り減らす（flag_quick=TRUE）
  • for文をベクトル演算で書ける部分はそうする（これやるとバグ取りが難しくなり得るので注意）。
• 実行時の動画。flag_quick=FALSEにした場合なので、TRUEにすれば更に２〜３倍の速さになると思う(ただしTRUEの場合は適応度関数の変更に対応していない)。
  
• プログラムの中身は以下の通り。

  K0=10.0; ##環境収容力の最大値
  sK=1.0; ##環境収容力の幅
  sa=0.5; ##競争効果の幅
  ##sa=0.7; 
  ##sa=1.2;
  ndiv=64; ##形質軸の刻み幅
  edge_die=1e-7; ##最小個体数密度。これ以下は絶滅。
  repeat_p=1000; ##変異体を出現させずに個体数動態をまわすステップ数（速く計算するため）
  dt=0.1; ##時間刻み幅
  mutation_rate=3.0e-4; ##突然変異率
  t_end=1500; ##総タイムステップ数
  set.seed(23); ##乱数の種。これをコメントアウトすれば毎回少し違う進化動態になる
  
  xx=seq(-1,1,length=ndiv);  ##とり得る形質値のセット
  n=xx*0.0;  ##形質軸上の個体数分布
  ancestor=c(ndiv/8); ##初期の表現型を指定。
  
  flag_quick=FALSE;
  
   ##環境収容力の分布
  K <- function(x){
      return (K0*exp(-x^2/(2.0*sK^2)));
  }
  
   ##競争効果の関数
  alpha <- function (x0,x1){
      return (exp(-(x0-x1)^2/(2.0*sa^2)));
  }
  
   ##変異型x1の適応度
  fitness <-function(x1,xx,n){
      return(1.0-sum(alpha(xx,x1)*n)/K(x1));
  }
  
   ##突然変異の関数
  mutate <- function(xx,n){
      for(i in 1:length(xx)){
          if(rbinom(1,size=1,prob=mutation_rate*n[i]*repeat_p*dt)>0){
              if(rbinom(1,size=1,prob=0.5)>0){
                  if((i < ndiv))n[i+1]=n[i+1]+edge_die*2.0;
              }else{
                  if((i > 1))n[i-1]=n[i-1]+edge_die*2.0;
              }
          }
      }
      return(n);
  }
  
  n[ancestor]=K(xx[ancestor])*0.5;  ##ancestorの位置に０でない個体数をセット
  n_out=matrix(ndiv*(t_end+1),nrow=ndiv,ncol=(t_end+1)); ##出力用行列（個体数分布の時系列データ）
  t_out=0.0*(1:(t_end+1)); ##出力用の時刻データ
  time=0.0;  ##時刻
  n_out[,1]=n;  ##出力用行列に初期状態を格納
  
  t_exec_start=proc.time();
  ##ここからメインループ
  for(t in 1:t_end){
  
      n=mutate(xx,n);  ##突然変異   
      mask_n=(n>edge_die);  ##存在する（最小個体数密度を上回る）表現型に目印をつける
      xxb=xx[mask_n];  ##存在する表現型だけにする（計算を速くするため）
      nb=n[mask_n];  ##存在する表現型だけの個体数密度
      dnb=nb*0.0;
      fit=dnb*0.0;
      if(flag_quick==TRUE){
          len_nb=length(nb);
          alphab=matrix(len_nb*len_nb,nrow=len_nb,ncol=len_nb);
          for(i in 1:len_nb){
              for(j in 1:len_nb)alphab[i,j]=alpha(xxb[i],xxb[j]);
          }
          Kb=K(xxb);
      }
          
      for(k in 1:repeat_p){
  
          if(flag_quick==TRUE){
              fit=(1.0-colSums(alphab*nb)/Kb); ##適応度の計算
              nb=nb+dt*nb*fit;  ##個体数密度を変化させる
          }
          else{
              for(i in 1:length(fit))fit[i]=fitness(xxb[i],xxb,nb);  ##適応度の計算 
              dnb=nb*fit;  ##増加率の計算
              nb=nb+dt*dnb; ##個体数密度を変化させる
          }
          
          time=time+dt;  ##時刻を更新
      }
      n[mask_n]=nb;  ##全ての表現型の個体数のベクトルに戻す
      n[which(n<=edge_die)]=0.0;  ##edge_die以下の表現型は絶滅
       
      n_out[,t+1]=n;  ##出力用行列に格納
      t_out[t+1]=time;  ##出力用の時刻配列に格納
  
      ##edge_die以上の表現型を1としてプロンプトに表示
      if(t%%5==1){
          cat("x: ",paste0(as.integer(n>edge_die),collapse=""));  
          cat("  time:", time,"\n");  
      }
  }
  
  cat("calculation time:\n");
  print(proc.time()-t_exec_start);
  ##10列おきにデータを使ってプロット
  list=which(seq(length(t_out))%%10==1);  
  filled.contour(xx,t_out[list],log(n_out[,list]+edge_die),xlab="Trait x",ylab="Time",main=paste0("sa:",sa, ", sK:",sK));
  ##dev.print(pdf,"test.pdf");  ##このコメントアウトを外せばtest.pdfとして画像が保存される（linuxだけ？）。
  

© 2020- Hiroshi C. Ito