clear all;
n_trial=10;%set the number of trials
%WT Cic intensity
load Cic_wt_dynamics.mat;
profile_exp_norm=9*profile_exp/(max(profile_exp));% Normalizing the WT Cic intensity
profile_exp_std=9*profile_exp_std/(max(profile_exp));% Normalizing the WT Cic intensity error
%Cic intensity in optogenetically activated embryo
load lighton_dynamics.mat;
profile_light_norm=c_n*profile_exp_norm(70)/c_n(1);%Normalizing the Cic intensity in optogenetically activated experiment with the wild type intensity at the same time
std_cic=std_cic*profile_exp_norm(70)/c_n(1);%Normalizing the error in the Cic intensity in optogenetically activated experiment
%Fraction of Cic bound in response to a step and pulse signal 
frac=[0.836957 0.099325 0.03036 0.0373 0.5641 1.1841];%Experimental estimation of fraction of Cic bound normalized such that the WT fraction at cycle 14 is 1.
frac_error=[0.1717 0.0819 0.0155 0.006 0.221 0.449];%Experimental error in estimation of fraction of Cic bound normalized such that the WT fraction at cycle 14 is 1.
% Transcription of tll  in response to a step and pulse signal 
load transcript_expt.mat;
trans_expt=cat (1,on_dynamics,on_off_dynamics);
trans_expt_error=cat (1,on_std,on_off_std);
% 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;%nuclear volume
     a_n(cycle)=pi*(d_n(cycle))^2;%nuclear surface area
     v_c(cycle)=(25 *sum(area)*(17.5*4.62)^2)/ncell-v_n(cycle);%cytoplasmic volume=average cellular volume-nuclear volume
     ncell=ncell*2;
end
trial=0;
while (trial<n_trial)
 
%initial guesses
k_im_0=10^(log10(0.1) + (log10(10)-log10(0.1))*rand);
k_ex_0=10^(log10(0.1) + (log10(10)-log10(0.1))*rand);
t_st=5+10*rand();
k_im_p=10^(log10(0.01) + (log10(1)-log10(0.01 ))*rand);
k_ex_p=10^(log10(0.01) + (log10(1)-log10(0.01))*rand);
k_p=(10^(log10(1) + (log10(50)-log10(1))*rand))/k_deg;
k_dp=(10^(log10(1) + (log10(50)-log10(1))*rand))/k_deg;
k_b=(10^(log10(0.01) + (log10(1)-log10(0.01))*rand))/k_deg;
k_ub=(10^(log10(0.01) + (log10(1)-log10(0.01))*rand))/k_deg;
k_on=10^(log10(0.01) + (log10(1)-log10(0.01))*rand)/k_deg;
k_off=10^(log10(0.01) + (log10(1)-log10(0.01))*rand)/k_deg;
m=20*rand();
%Find error at initial guesses
    simulate_unphoshpho;%Simulation for nucleus at the middle of the embryo
    error_1=norm (profile_simu_norm-profile_exp_norm)/norm (profile_exp_std);
    
    pole_simulate;%Simulation for a nucleus at the pole
    profile_pole=profile_pole/profile_c (round ((25+t_st)*2),1);
    error_2=abs (2.8-profile_simu_norm(round(50+2*t_st),1)/profile_pole (round (t_st*2),1))/0.15;
    
    simulate_light_phoshpho;%Simulation for nucleus at the middle of the embryo optogenetically activated at cycle 14  
    simulate_light_on_off;%Simulation for nucleus at the middle of the embryo optogenetically activated at cycle 14  but turned off after 5 mins
    profile_simu_light_norm=profile_simu_light/profile_c (round ((25+t_st)*2),1);
    error_3=norm (profile_simu_light_norm(71:129)-profile_light_norm)/norm (std_cic); %error from the dynamics of Cic optogenetically active embryo
    frac_simu=[profile_b_on(22,1) profile_b_on(30,1) profile_b_on(82,1) profile_b_on_off(40,1) profile_b_on_off(50,1) profile_b_on_off(60,1)];
    error_4=norm (frac_simu-frac)/norm (frac_error);%error from the immunoprecipiation data
    transcript_simu=cat (1,profile_t_on(20:50),profile_t_on_off(10:50));
    error_5=norm (transcript_simu-trans_expt)/norm (trans_expt_error);%error from the transcriptional response

    prev_error=[error_1 error_2 error_3 error_4 error_5];
    prev_para=[k_im_0 k_ex_0 t_st k_im_p k_ex_p k_p k_dp k_b*c_max k_ub k_on k_off m];
%monte carlo parameters
Delta=0.1;
Beta=500;


for (step=1:2000)
    v=randn(12,1);
    change=zeros(12);
    for (i=1:12) 
        change(i,i) =1+Delta*v(i,1)/norm(v); 
    end
    new_para=prev_para*change;
    %updated parameters
    k_im_0=new_para(1); k_ex_0=new_para(2); t_st=new_para(3); 
    k_im_p=new_para(4); k_ex_p=new_para(5); 
    k_p=new_para(6); k_dp=new_para(7); 
    k_b=new_para(8); k_ub=new_para(9);
    k_on=new_para(10); k_off=new_para(11); m=new_para(12); 
    
    simulate_unphoshpho;%Simulation for nucleus at the middle of the embryo
    error_1=norm (profile_simu_norm-profile_exp_norm)/norm (profile_exp_std);
    
    pole_simulate;%Simulation for a nucleus at the pole
    profile_pole=profile_pole/profile_c (round ((25+t_st)*2),1);
    error_2=abs (2.8-profile_simu_norm(round(50+2*t_st),1)/profile_pole (round (t_st*2),1))/0.15;
    
    simulate_light_phoshpho;%Simulation for nucleus at the middle of the embryo optogenetically activated at cycle 14  
    simulate_light_on_off;%Simulation for nucleus at the middle of the embryo optogenetically activated at cycle 14  but turned off after 5 mins
    profile_simu_light_norm=profile_simu_light/profile_c (round ((25+t_st)*2),1);
    error_3=norm (profile_simu_light_norm(71:129)-profile_light_norm)/norm (std_cic); %error from the dynamics of Cic optogenetically active embryo
    frac_simu=[profile_b_on(22,1) profile_b_on(30,1) profile_b_on(82,1) profile_b_on_off(40,1) profile_b_on_off(50,1) profile_b_on_off(60,1)];
    error_4=norm (frac_simu-frac)/norm (frac_error);%error from the immunoprecipiation data
    transcript_simu=cat (1,profile_t_on(20:50),profile_t_on_off(10:50));
    error_5=norm (transcript_simu-trans_expt)/norm (trans_expt_error);%error from the transcriptional response
    
    new_error=[error_1 error_2 error_3 error_4 error_5];
    accept=0;
        if (norm(new_error)< norm(prev_error))
            accept=1;%accept the move if the error reduces

        else
            prob=exp(-Beta*((norm(new_error))^2-(norm(prev_error))^2));
            if (rand()<prob) accept=1;end%accept the move with a probability P if the error increases
        end
        if (accept==1)%if accepted update the parameters
          prev_error=new_error;
          prev_para=new_para;
        end
   
end
  
  if (max (prev_error) >1.5) continue; end;%Reject the trials resulting in huge errors
  trial=trial+1;
  para_trial_final (trial,1:12)=prev_para;
  error_trial_final (trial,1:5)=prev_error;
  
  save ('MCMC_para_demo.mat','para_trial_final','error_trial_final');%save the parameters
end
