clear all;
% known parameters
k_deg=1.4/15;%rate of degradation of Cic
k_l_on=(1/(0.01*k_deg));%rate of optogenetic activation of ERK
k_l_off=(1/(2*k_deg));%rate of deactivation of ERK
E_max=2;%The total pool of ERK in comparison to the ERK level at the pole of the WT
cycletime=[5 7.5 12.5 45]*k_deg;%Duration of nuclear cycle
d_n=[10.5 9.2 8.2 6.5];%Diameter of nuclei
ncell=750;%initial number of nuclei at cycle 11.

for (cycle=1:4)
     load (strcat('voro_1',num2str(cycle),'.mat'));
      v_n(cycle)=pi*((d_n(cycle))^3)/6;
      a_n(cycle)=pi*(d_n(cycle))^2;
      v_c(cycle)=(25 *sum(area)*(17.5*4.62)^2)/ncell-v_n(cycle);
     ncell=ncell*2;
end
load MCMC_para.mat;
for (trial_n=1:size(para_trial_final,1))
    %load the unknown parameters
    k_im_0=para_trial_final (trial_n,1);
    k_ex_0=para_trial_final (trial_n,2);
    t_st=para_trial_final (trial_n,3);
    k_im_p=para_trial_final (trial_n,4);
    k_ex_p=para_trial_final (trial_n,5);
    k_p=para_trial_final (trial_n,6);
    k_dp=para_trial_final (trial_n,7);
    k_b= para_trial_final (trial_n,8);
    k_ub=para_trial_final (trial_n,9);
    k_on=para_trial_final (trial_n,10);
    k_off=para_trial_final (trial_n,11);
    m=para_trial_final (trial_n,12);
    % Nuclei at the middle of the wild type embryo 
    simulate_unphoshpho;
    nuclear_wt (trial_n,:)=profile_simu_norm;
    % Nuclei at the middle of an wild type embryo activated with a step signal
    simulate_light_phoshpho; 
    nuclear_on(trial_n,:)=profile_simu_light/profile_c (round ((25+t_st)*2),1);
    bound_on(trial_n,:)=100*profile_b_on(21:60);
    transcript_on(trial_n,:)=100*profile_t_on(21:60);
    % Nuclei at the middle of an wild type embryo activated with a pulse signal
    simulate_light_on_off;
    bound_on_off(trial_n,:)=100*profile_b_on_off(21:60);
    transcript_on_off(trial_n,:)=100*profile_t_on_off(21:60);
    % Nuclei at the pole of an wild type embryo 
    pole_simulate;
    profile_pole=profile_pole/profile_c (round ((25+t_st)*2),1);
    middle_pole_ratio(trial_n)=profile_simu_norm(round(50+2*t_st),1)/profile_pole (round (t_st*2),1);
    save ('ensemble_data.mat','nuclear_wt','nuclear_on','bound_on','transcript_on','bound_on_off','transcript_on_off','middle_pole_ratio');
    disp (trial_n);
end

