How to implement mass balance in ODE solver?

I am modelling the chemical kinetics of a radiolysis system. My code works fine but the results do not match the timescale of my experimental data, and I am wondering if this is because my current implementation of the ODE solver does not handle mass balance in a way that reflects conservation of mass in a physical system. If this is correct, is there a way to implement mass balance in the ODE solver?

To further explain my question: in the physical system I am modelling, the rate of concentration change for all species is limited by the concentration of one species, the hydroxyl radical OH, which is generated by water radiolysis at a rate of 1E-6 M/s. At t=0, the only species present in the system are OH and acetate (initial conc 1E-3 M), and the model fits the pseudo-first order curve seen in the experimental data. As the system evolves with time, more organic species are formed which can also react with OH. Because OH is at a very low steady state concentration, there is competition between these species and acetate for reaction with OH, thereby reducing the rate of acetate decomposition. However, as rate of acetate consumption changes very little in my model, I believe the model is not taking into account this competition, and instead behaving as though the full concentration of OH at time t is available to every species in the system, therefore violating the principle of conservation of mass.

I have included a shortened version of my current model here which should produce a graph of results; in the model, acetate is fully consumed around 1400 s, whereas in experimental data this is closer to 10,000 s. I have a symbolic expression of the ODE system as well if that is helpful.

A=[-1	1	0	0	0	-1	0	0	0	0	0
0	-2	0	0	0	0	0	1	0	0	0
0	0	0	0	0	-1	0	-1	1	0	0
0	0	0	1	0	-1	0	0	-1	0	0
0	-1	0	0	0	0	0	0	-1	1	0
0	0	0	0	0	-1	0	0	0	-1	1
0	0	1	0	0	-1	0	0	0	0	-1
0	0	0	0	0	1	0	0	0	0	0
0	0	0	-1	0	-1	1	0	0	0	0
0	0	0	0	1	-1	-1	0	0	0	0];
k=[85000000
520000000
310000000
860000000
10000000
100000000
800000000
1.46000000000000e-06
10000000
10000000];
species=["Acetate","AcetateR","Citrate","Malate","Malonate","OH","Oxaloacetate","Succinate","SuccinateR","Tricarballylate","TricarballylateR"];
options=odeset('NonNegative',1,'RelTol',1E-14,'AbsTol',1E-17);
cinit=zeros(size(A,2),1); cinit(1)=1E-3;
tspan=[0 1434];
[t,Y]=ode15s(@(t,conc) BUILD_RATES(k,A,conc,t),tspan,cinit,options);

% SEPARATE RADICALS FROM NON RADICALS, PLOT NON RADICAL SPECIES
[YnonRad,nonRadicals,radicals]=SEPARATE_RADICALS(Y,species);
plot(t,YnonRad); legend(nonRadicals)

%% Build the numerical solver for the ODE system, returns dYdt for use with ode15s, requires k, stoichmatrix inputs
function [dYdt]=BUILD_RATES(k,stoichmatrix,conc,~)
ratevector=zeros(length(k),1); sz=size(stoichmatrix); nrxns=sz(1); n=sz(2); 
for i=1:nrxns % Loop over all reactions
    ratevector(i)=k(i); % Begin each line of rate eqn vector with k value for that reaction
    for j=1:n % Loop over all species
        if stoichmatrix(i,j)<0 % Reactants have negative values on the stoichiometric matrix
            ratevector(i)=ratevector(i)*(conc(j)^(abs(stoichmatrix(i,j)))); 
            % Multiply current rate vector by next reactant
            % (to the power of the positive of its reaction order)
        end
    end
end
dYdt=stoichmatrix'*ratevector; % Create the concentration change matrix by  
                               % multiplying rate vector by stoichmatrix transpose
end
%% Separate radical and non-radical species to allow for clearer plots
function [Ynonrad,nonRadicals,radicals]=SEPARATE_RADICALS(Y,species)
j=0; k=0;
for i=1:length(species)
    if regexp(species(i),'.+R')
        j=j+1;
        radicals(j)=species(i);
    else
        k=k+1;
        nonRadicals(k)=species(i);
        Ynonrad(:,k)=Y(:,i);
    end
end
end

[Bjorn Gustavsson]

Can you also write the chemical reactions or the species continuity-equations? That would help reading your code since that way we can directly see what reactions are losses and sources of what substances. Not everyone has the chemical konwledge to find problems with imballances of production/losses of for example "Oxaloacetate"...

Answer 1

[Steven Lord]

I think you could adapt the Robertson problem example from this documentation page.