 c_u=zeros (8,14000,4);
 c_p=zeros (8,14000,4);
 c_n_u=zeros (8,14000,4);
 c_n_p=zeros (8,14000,4);
 B_u=zeros (8,14000,4);
 B_p=zeros (8,14000,4);
 f_t=zeros (8,14000,4);
  ncell=1;

                                for (cycle=1:4)

                                    if (cycle~=1)
                                        for (i=1:ncell/2)
                                            c_u(i,1,cycle)=c_next_u(i,1)/(2*v_c(cycle));
                                            c_u(0.5*ncell+i,1,cycle)=c_u(i,1,cycle);
                                            c_p(i,1,cycle)=c_next_p(i,1)/(2*v_c(cycle));
                                            c_p(0.5*ncell+i,1,cycle)=c_p(i,1,cycle);

                                        end
                                    end
                                for (i=1:ncell)
                                    c_n_u(i,1,cycle)=0;    
                                    c_n_p(i,1,cycle)=0;  
                                end
                                
                                dt=0.001;
                                time=cycletime (cycle)/dt;
                                count=0;
                                 k_in_p =k_im_p*(a_n(cycle)/a_n(3))*(v_n(3)/v_c(cycle))/k_deg;
                                 k_out_p=k_ex_p*(a_n(cycle)/a_n(3))*(v_n(3)/v_c(cycle))/k_deg; 
                                 k_in_u =k_im_0*(a_n(cycle)/a_n(3))*(v_n(3)/v_c(cycle))/k_deg;
                                 k_out_u =k_ex_0*(a_n(cycle)/a_n(3))*(v_n(3)/v_c(cycle))/k_deg;  
                                t_on=10;
                                 
                                 
                                for (t=1:time)
                                    for (i=1:ncell)
                                        if (cycle==4 && t*dt/k_deg >t_st) q=0;
                                        else q=1;
                                        end
                                        
                                        if (cycle==4 && t*dt/k_deg >t_on) 
                                            dperk(i,1)=dperk(i,1)+k_l_on*(E_max-dperk(i,1))*dt-k_l_off*dperk(i,1)*dt; 
                                        else
                                            dperk(i,1)=0;
                                        end
                                        c_u(i,t+1,cycle)=c_u(i,t,cycle)+(q-c_u(i,t,cycle)+k_out_u*c_n_u(i,t,cycle)-k_in_u*c_u(i,t,cycle)-k_p*dperk(i,1)*c_u(i,t,cycle)+k_dp*c_p(i,t,cycle))*dt;
                                        c_p(i,t+1,cycle)=c_p(i,t,cycle)+(-c_p(i,t,cycle)+k_out_p*c_n_p(i,t,cycle)-k_in_p*c_p(i,t,cycle)+k_p*dperk(i,1)*c_u(i,t,cycle)-k_dp*c_p(i,t,cycle))*dt;
                                        c_n_u(i,t+1,cycle)=c_n_u(i,t,cycle)+((-k_out_u*c_n_u(i,t,cycle)+k_in_u*c_u(i,t,cycle))*(v_c(cycle)/v_n(cycle))-k_p*dperk(i,1)*c_n_u(i,t,cycle)+k_dp*c_n_p(i,t,cycle)-c_max*0.1*(k_b*(1-B_u(i,t,cycle)-B_p(i,t,cycle))*c_n_u(i,t,cycle)-k_ub*B_u(i,t,cycle)))*dt;
                                        c_n_p(i,t+1,cycle)=c_n_p(i,t,cycle)+((-k_out_p*c_n_p(i,t,cycle)+k_in_p*c_p(i,t,cycle))*(v_c(cycle)/v_n(cycle))+k_p*dperk(i,1)*c_n_u(i,t,cycle)-k_dp*c_n_p(i,t,cycle)-c_max*0.1*((k_b/1000)*(1-B_u(i,t,cycle)-B_p(i,t,cycle))*c_n_p(i,t,cycle)-k_ub*B_p(i,t,cycle)))*dt;
                                        B_u(i,t+1,cycle)= B_u(i,t,cycle)+(-k_p*dperk(i,1)*B_u(i,t,cycle)+k_dp*B_p(i,t,cycle)+k_b*(1-B_u(i,t,cycle)-B_p(i,t,cycle))*c_n_u(i,t,cycle)-k_ub*B_u(i,t,cycle))*dt;
                                        B_p(i,t+1,cycle)= B_p(i,t,cycle)+(k_p*dperk(i,1)*B_u(i,t,cycle)-k_dp*B_p(i,t,cycle)+(k_b/1000)*(1-B_u(i,t,cycle)-B_p(i,t,cycle))*c_n_p(i,t,cycle)-k_ub*B_p(i,t,cycle))*dt;
                                       f_t(i,t+1,cycle)=f_t(i,t,cycle)+(k_on*(1-B_u(i,t,cycle)-B_p(i,t,cycle))^m*(1-f_t(i,t,cycle))-k_off*f_t(i,t,cycle))*dt; 
                                    end


                                    for (i=1:ncell)
                                        c_next_u(i,1)=c_n_u(i,t+1,cycle)*v_n(cycle) + c_u(i,t+1,cycle)* v_c(cycle);
                                        c_next_p(i,1)=(c_n_p(i,t+1,cycle))*v_n(cycle) + c_p(i,t+1,cycle)* v_c(cycle);
                                    end
                                end
                                 ncell=ncell*2;

end
                                
                                
                                profile_simu_light(1:10,1)=mean (c_n_u(1:1,round ((0.5:0.5:5)*k_deg/dt),1)+c_n_p(1:1,round ((0.5:0.5:5)*k_deg/dt),1)+c_max*0.1*(B_u(1:1,round ((0.5:0.5:5)*k_deg/dt),1)+B_p(1:1,round ((0.5:0.5:5)*k_deg/dt),1)),1);
                                profile_simu_light(11:25,1)=mean (c_n_u(1:2,round ((0.5:0.5:7.5)*k_deg/dt),2)+c_n_p(1:2,round ((0.5:0.5:7.5)*k_deg/dt),2)+c_max*0.1*(B_u(1:2,round ((0.5:0.5:7.5)*k_deg/dt),2)+B_p(1:2,round ((0.5:0.5:7.5)*k_deg/dt),2)),1);
                                profile_simu_light(26:50,1)=mean (c_n_u(1:4,round ((0.5:0.5:12.5)*k_deg/dt),3)+c_n_p(1:4,round ((0.5:0.5:12.5)*k_deg/dt),3)+c_max*0.1*(B_u(1:4,round ((0.5:0.5:12.5)*k_deg/dt),3)+B_p(1:4,round ((0.5:0.5:12.5)*k_deg/dt),3)),1);
                                profile_simu_light (51:140,1)=mean (c_n_u(1:8,round ((0.5:0.5:45)*k_deg/dt),4)+c_n_p(1:8,round ((0.5:0.5:45)*k_deg/dt),4)+c_max*0.1*(B_u(1:8,round ((0.5:0.5:45)*k_deg/dt),4)+B_p(1:8,round ((0.5:0.5:45)*k_deg/dt),4)),1);
                                profile_b_on (1:90,1)=mean ((B_u(1:8,round ((0.5:0.5:45)*k_deg/dt),4)+B_p(1:8,round ((0.5:0.5:45)*k_deg/dt),4)),1);
                                profile_t_on (1:90,1)=mean (f_t(1:8,round ((0.5:0.5:45)*k_deg/dt),4),1);
                     