function [X,t] = BK_stochastic(X0,t_end,dt,struct)
%% Stochastic opening of BK channel (10 states)
%Dimension of the state space
dimX = size(X0,1);
t = 0:dt:t_end;% ms
allocated_memory = length(t);

%initialization of the states
MC0 = zeros(1,allocated_memory);
MC1 = zeros(1,allocated_memory);
MC2 = zeros(1,allocated_memory);
MC3 = zeros(1,allocated_memory);
MC4 = zeros(1,allocated_memory);
MC5 = zeros(1,allocated_memory);
MC6 = zeros(1,allocated_memory);
MC7 = zeros(1,allocated_memory);
MC8 = zeros(1,allocated_memory);
MC9 = zeros(1,allocated_memory);
MC0(1) = X0(1);
MC1(1) = X0(2);
MC2(1) = X0(3);
MC3(1) = X0(4);
MC4(1) = X0(5);
MC5(1) = X0(6);
MC6(1) = X0(7);
MC7(1) = X0(8);
MC8(1) = X0(9);
MC9(1) = X0(10);
%% reactions
s_to_ms = 1e-3;
M_to_microM = 1e-6;


KonO = 1e9*s_to_ms*M_to_microM;
KonC = 1e9*s_to_ms*M_to_microM;
KoffO = 1065*s_to_ms;
KoffC = 11917*s_to_ms;


nReactions=26;
Caloc = struct.Ca;
%rng(seed_number)
for ind=1:allocated_memory
  % random variables
    V = interp1(struct.Voltage_time,struct.Voltage,t(ind));
    Y = rand(1);
    Q = zeros(1,nReactions);
    P0 = MC0(ind);
    P1 = MC1(ind);
    P2 = MC2(ind);
    P3 = MC3(ind);
    P4 = MC4(ind);
    P5 = MC5(ind);
    P6 = MC6(ind);
    P7 = MC7(ind);
    P8 = MC8(ind);
    P9 = MC9(ind);
    
    %% time dependent variables
    
F = 96485.34;
R = 8.314;
T = 310; 
qf = 0.73;
qb = 0.58;
esp=(F*V*1e-3)/(R*T);

alfa1 = 5.5*s_to_ms*exp(qf*esp);
beta1 = 8669*s_to_ms*exp(-qb*esp);
alfa2 = 8*s_to_ms*exp(qf*esp);
beta2 = 1127*s_to_ms*exp(-qb*esp);
alfa3 = 2*s_to_ms*exp(qf*esp);
beta3 = 25.2*s_to_ms*exp(-qb*esp);
alfa4 = 884*s_to_ms*exp(qf*esp);
beta4 = 1013*s_to_ms*exp(-qb*esp);
alfa5 = 900*s_to_ms*exp(qf*esp);
beta5 = 125.7*s_to_ms*exp(-qb*esp);

    
    
    KO=1.06;
    KC=11.92;
    
    %% Propensity of the controller----------------------------------------
    Q(1,1) = 4*KonO*Caloc(ind)*dt*P0;
    Q(1,2) = beta1*dt*P0; 
    Q(1,3) = 3*KonO*Caloc(ind)*dt*P1;
    Q(1,4) = beta2*dt*P1; 
    Q(1,5) = KoffO*dt*P1;
    Q(1,6) = 2*KonO*Caloc(ind)*dt*P2;
    Q(1,7) = beta3*dt*P2;
    Q(1,8) = 2*KoffO*dt*P2;
    Q(1,9) = KonO*Caloc(ind)*dt*P3;
    Q(1,10) = beta4*dt*P3;
    Q(1,11) = 3*KoffO*dt*P3;
    Q(1,12) = beta5*dt*P4;
    Q(1,13) = 4*KoffO*dt*P4;
    Q(1,14) = 4*KonC*Caloc(ind)*dt*P5;
    Q(1,15) = alfa1*dt*P5;
    Q(1,16) = 3*KonC*Caloc(ind)*dt*P6;
    Q(1,17) = alfa2*dt*P6;
    Q(1,18) = KoffC*dt*P6;
    Q(1,19) = 2*KonC*Caloc(ind)*dt*P7;
    Q(1,20) = alfa3*dt*P7;
    Q(1,21) = 2*KoffC*dt*P7;
    Q(1,22) = KonC*Caloc(ind)*dt*P8;
    Q(1,23) = alfa4*dt*P8;
    Q(1,24) = 3*KoffC*dt*P8;
    Q(1,25) = alfa5*dt*P9;
    Q(1,26) = 4*KoffC*dt*P9;
    P=[P0;P1;P2;P3;P4;P5;P6;P7;P8;P9];
   
    current_state = find(P==1);
    
    if current_state == 1 %MC0
    
        if Y < Q(1,1) %4*KonO*KO*P0
            MC0(ind+1) = MC0(ind)-1;
            MC1(ind+1) = MC1(ind)+1;
       elseif Y < Q(1,1)+Q(1,2)% alfa1*P0
            MC0(ind+1) = MC0(ind)-1;
            MC5(ind+1) = MC5(ind)+1;
        else 
            MC0(ind+1) = MC0(ind);   
        end
        
    elseif current_state == 2 %MC1
        
        if Y < Q(1,3)%3*KonO*KO*P1
            MC1(ind+1) = MC1(ind)-1;
            MC2(ind+1) = MC2(ind)+1;
       elseif Y < Q(1,3)+Q(1,4) %beta1*P1
            MC1(ind+1) = MC1(ind)-1;
            MC6(ind+1) = MC6(ind)+1;
       elseif Y < Q(1,3)+Q(1,4)+Q(1,5) %KoffO*P1
            MC0(ind+1) = MC0(ind)+1;
            MC1(ind+1) = MC1(ind)-1;
        else
            MC1(ind+1) = MC1(ind);
        end
   
     elseif current_state == 3 %MC2
         
         if Y<Q(1,6)% 2*KonO*KO*P2
            MC2(ind+1) = MC2(ind)-1;
            MC3(ind+1) = MC3(ind)+1;
        elseif Y < Q(1,6)+Q(1,7) %gamma1*P2
            MC2(ind+1) = MC2(ind)-1;
            MC7(ind+1) = MC7(ind)+1;
         elseif Y < Q(1,6)+Q(1,7)+Q(1,8) %2*KoffO*P2
            MC2(ind+1) = MC2(ind)-1;
            MC1(ind+1) = MC1(ind)+1;
         else
            MC2(ind+1) = MC2(ind);
         end
         
     elseif current_state == 4 %MC3
         
         if Y < Q(1,9) %KonO*KO*P3
            MC4(ind+1) = MC4(ind)+1;
            MC3(ind+1) = MC3(ind)-1;
        elseif Y < Q(1,9)+Q(1,10) %delta1*P3
            MC3(ind+1) = MC3(ind)-1;
            MC8(ind+1) = MC8(ind)+1;
         elseif Y < Q(1,9)+Q(1,10)+Q(1,11) % 3*KoffO*P3
            MC3(ind+1) = MC3(ind)-1;
            MC2(ind+1) = MC2(ind)+1;
         else
            MC3(ind+1) = MC3(ind);
         end
        
     elseif current_state == 5 %MC4
         
         if Y<Q(1,12) %csi1*P4
             MC9(ind+1) = MC9(ind)+1;
             MC4(ind+1) = MC4(ind)-1;
        elseif Y < Q(1,12)+Q(1,13) % 4*KoffO*P4
             MC4(ind+1) = MC4(ind)-1;
             MC3(ind+1) = MC3(ind)+1;
         else
             MC4(ind+1) = MC4(ind);
         end
        
          elseif current_state == 6 %MC5
              
        if Y < Q(1,14)  %4*KonC*KC*P5
             MC6(ind+1) = MC6(ind)+1;
             MC5(ind+1) = MC5(ind)-1;
        elseif Y < Q(1,14)+Q(1,15) %alfa2*P5
             MC5(ind+1) = MC5(ind)-1;
             MC0(ind+1) = MC0(ind)+1;
        else
             MC5(ind+1) = MC5(ind);
        end 
              
        elseif current_state == 7
           if Y < Q(1,16) %3*KonC*KC*P6
             MC7(ind+1) = MC7(ind)+1;
             MC6(ind+1) = MC6(ind)-1;
           elseif Y < Q(1,16)+Q(1,17) %beta2*P6
             MC6(ind+1) = MC6(ind)-1;
             MC1(ind+1) = MC1(ind)+1;
           elseif Y < Q(1,16)+Q(1,17)+Q(1,18) %KoffC*P6
             MC6(ind+1) = MC6(ind)-1;
             MC5(ind+1) = MC5(ind)+1;
           else
             MC6(ind+1) = MC6(ind);
           end
           
           elseif current_state == 8 %MC7
               
           if Y < Q(1,19) %2*KonC*KC*P7
              MC8(ind+1) = MC8(ind)+1;
              MC7(ind+1) = MC7(ind)-1;
           elseif Y<Q(1,19)+Q(1,20) %gamma2*P7
              MC7(ind+1) = MC7(ind)-1;
              MC2(ind+1) = MC2(ind)+1;
           elseif Y<Q(1,19)+Q(1,20)+Q(1,21) %2*KoffC*P7
              MC7(ind+1) = MC7(ind)-1;
              MC6(ind+1) = MC6(ind)+1;
           else
              MC7(ind+1) = MC7(ind);     
           end
           
           elseif current_state==9 %MC8
               
           if Y<Q(1,22) %KonC*KC*P8
              MC9(ind+1) = MC9(ind)+1;
              MC8(ind+1) = MC8(ind)-1;
           elseif Y<Q(1,22)+Q(1,23) %delta2*P8
              MC8(ind+1) = MC8(ind)-1;
              MC3(ind+1) = MC3(ind)+1;
           elseif Y<Q(1,22)+Q(1,23)+Q(1,24) %3*KoffC*P8
              MC8(ind+1) = MC8(ind)-1;
              MC7(ind+1) = MC7(ind)+1;      
           else
              MC8(ind+1) = MC8(ind);   
           end

           elseif current_state==10 %MC9
               
           if Y<Q(1,25) %csi2*P9
              MC4(ind+1)=MC4(ind)+1;
              MC9(ind+1)=MC9(ind)-1;
           elseif Y<Q(1,25)+Q(1,26)%4*KoffC*P9
              MC9(ind+1)=MC9(ind)-1;
              MC8(ind+1)=MC8(ind)+1;
           else
              MC9(ind+1)=MC9(ind);
              
           end           
    end

     
end
X(1,:) = MC0(1:allocated_memory);
X(2,:) = MC1(1:allocated_memory);
X(3,:) = MC2(1:allocated_memory);
X(4,:) = MC3(1:allocated_memory);
X(5,:) = MC4(1:allocated_memory);
X(6,:) = MC5(1:allocated_memory);
X(7,:) = MC6(1:allocated_memory);
X(8,:) = MC7(1:allocated_memory);
X(9,:) = MC8(1:allocated_memory);
X(10,:) = MC9(1:allocated_memory);


