function dcdt = IsoEquilibration_JW(t, C)

global  H   OH   x1 
global  k1f   k1r   k4f   k4r
global  x1_13   x1_18   x1_47
global  k1f_13   k1r_13   k4f_13   k4r_13 
global  a1f   a1r   a4f   a4r
global  b1f   b1r   b4f   b4r 
global  R18w   R18OH
global  s1f   s1r   s4f   s4r
global  p1f   p1r   p4f   p4r

    co2d    = 1;        
    DCP     = 2;        
    co2d_13 = 3;        
    DCP_13  = 4;        
    co2d_18 = 5;       
    DCP_18  = 6;        
    co2d_47 = 7;
    DCP_47  = 8; 


dcdt(co2d) = -(k1f + k4f*OH) * C(co2d) ...
    +(k1r*H + k4r) * C(DCP) * x1;

dcdt(DCP) = (k1f + k4f*OH) * C(co2d) ...
    -(k1r*H + k4r) * C(DCP) * x1;

dcdt(co2d_13) = -(k1f_13 + k4f_13*OH)* C(co2d_13) ...
    +(k1r_13*H + k4r_13) * C(DCP_13) * x1_13; 

dcdt(DCP_13) = (k1f_13 + k4f_13*OH) * C(co2d_13) ...
    -(k1r_13*H + k4r_13) * C(DCP_13) * x1_13; 

dcdt(co2d_18) = -(b1f + b4f*OH) * C(co2d_18) ...
    + (2/3) * (b1r*H + b4r) * C(DCP_18) * x1_18;  
                                     % After Christensen et al. 2021

dcdt(DCP_18) = (b1f + b4f*OH) * C(co2d_18) ...
    -(2/3) * (b1r*H + b4r) * C(DCP_18) * x1_18 ...
    + (a1f*R18w + a4f*R18OH*OH) * C(co2d) ...
    -(1/3) * (a1r*H + a4r) * C(DCP_18) * x1_18;   
                                     % After Christensen et al. 2021

dcdt(co2d_47) = -(s1f + s4f*OH) * C(co2d_47) ...
    + (2/3) * (s1r*H + s4r) * C(DCP_47) * x1_47;
                                    % following Christensen et al. 2021
                                    % for oxygen

dcdt(DCP_47) = (s1f + s4f*OH) * C(co2d_47) ...
    -(2/3) * (s1r*H + s4r) * C(DCP_47) * x1_47 ...
    + (p1f*R18w + p4f*R18OH*OH) * C(co2d_13) ...
    -(1/3) * (p1r*H + p4r) * C(DCP_47) * x1_47;
                                    % following Christensen et al. 2021
                                    % for oxygen

dcdt = (dcdt)';