How to use ode15s with a constant that varies by two parameters?

I am using Matlab to simulate protein electrokinetics. Part of the complexity of simulating a system is that some of the constants vary by a time independent d parameter and a time dependent parameter. I have been able to incorporate the time dependency of the constant correctly, however, the time independent parameter is stumping me. I tried to define the anonymous function I used for the time dependent part as having two variables and then giving it the values I want the d parameter to go through and then taking the average, but this is not giving me the desired results. The code:

Enfcn = @(t) Ei + v * t;
kred = @(t, d) 225000 * exp(-d) * exp(((-a * F)/(R * T)) * (Enfcn(t) - Eo));
kox  = @(t, d) 225000 * exp(-d)* exp((((1-a) * F)/(R * T)) * (Enfcn(t) - Eo));

C0 = [0.01, 0, 0, 0, 0, 0, 1e9, 0];

options = odeset('RelTol', 1e-9, 'AbsTol', 1e-10);

[Time, Concentrations] = ode15s(@(t, C) mechanism(t, C, kred, kox), tspan, C0, options);

function dC = mechanism(t, C, kred, kox)

 x = 4:0.1:12;
     % Rate Laws
        dNiII   = -y(1)*y(7) + y(3) + y(2)*y(7);
        dNiIISH = -y(2) + y(2)*y(7) - kred(t, x)*y(2) + kox(t, x)*y(3);
        dNiISH  = kred(t, x)*y(2) - kox(t, x)*y(3) - y(3);
        de      = -kred(t, x)*(y(2) + y(5)) + kox(t, x)*(y(3) + y(6));
        dNiIIIH = y(3) - kred(t, x)*y(5) + kox(t, x)*y(6);
        dNiIIH  = kred(t, x)*y(5) - kox(t, x)*y(6) - y(7)*y(8);
        dH      = -y(6)*y(7) - y(1)*y(7) + kb*y(6);
        dH2     = y(6)*y(7);

     % Assign Output Variables
        dC(1,:) = mean(dNiII);
        dC(2,:) = mean(dNiIISH);
        dC(3,:) = mean(dNiISH);
        dC(4,:) = mean(dNiIIIH);
        dC(5,:) = mean(dNiIIH);
        dC(6,:) = mean(de)*(96485.33289/1e9);
        dC(7,:) = mean(dH);
        dC(8,:) = mean(dH2);
end

I am intending for the differential equations to be solved for time, and then be solved as a function of the d parameter. I was thinking of using nested functions but I am not sure how to execute that. At any rate, for a particular time, the constants kred and kox should have say 10 values based on the d parameters. I want for the system to be solved with one constant of d for time and then all the solutions for all the different d values are averaged, but the above code is not accomplishing what I hoped. Thanks for any help in advance.

Answer 1

[J. Alex Lee]

From you description, it sounds like you want to do the average over d-values outside of our ODE solution. But in your code, you are averaging for your d-values inside the odefun, which means you are averaging your differential equations, not your solutions.

If this is the case, define the odefun for a single value of d, and do your loop outside it:

% define list of d values 
x = 4:0.1:12;

% loop over desired list of d values
for i = 1:length(x)
    d = x(i)
    [Time{i}, Concentrations{i}] = ode15s(@(t, C) mechanism(t, C, kred, kox,d), tspan, C0, options);
end

M = length(x)
AvgConc = zeros(length(tspan),length(C0));
% average your solutions
for i = 1:M
    AvgConc = AvgConc + Concentrations{i}/M
end

% odefun definition
function dC = mechanism(t, C, kred, kox , x) % <=== add d to the parameter ; i call it x because you used x below
     % Rate Laws
        dNiII   = -y(1)*y(7) + y(3) + y(2)*y(7);
        dNiIISH = -y(2) + y(2)*y(7) - kred(t, x)*y(2) + kox(t, x)*y(3);
        dNiISH  = kred(t, x)*y(2) - kox(t, x)*y(3) - y(3);
        de      = -kred(t, x)*(y(2) + y(5)) + kox(t, x)*(y(3) + y(6));
        dNiIIIH = y(3) - kred(t, x)*y(5) + kox(t, x)*y(6);
        dNiIIH  = kred(t, x)*y(5) - kox(t, x)*y(6) - y(7)*y(8);
        dH      = -y(6)*y(7) - y(1)*y(7) + kb*y(6);
        dH2     = y(6)*y(7);
     % Assign Output Variables
        dC(1,:) = (dNiII);
        dC(2,:) = (dNiIISH);
        dC(3,:) = (dNiISH);
        dC(4,:) = (dNiIIIH);
        dC(5,:) = (dNiIIH);
        dC(6,:) = (de)*(96485.33289/1e9);
        dC(7,:) = (dH);
        dC(8,:) = (dH2);
end

[Riley Stein]

This is giving me answers closer to what's expected, but not exactly the values (magnitude off). A problem I have with both yours and what I tried is that the plots are too "linear" whereas what I'm trying to simulate bows out around a specific time and then levels off.

[J. Alex Lee]

That sounds like an application/physics/math problem rather than a Matlab/coding problem...

I don't think you're likely to get much more help from this community if that's the case.

[Riley Stein]

I see. I have a Mathematica script that produces expected values. I am new to Matlab and so I am having trouble transferring everything and I'm not sure what's going wrong where. Thank you for the help though. It's definitely closer to what I am expecting.

[J. Alex Lee]

if your equations are confirmed to be transcribed correctly from the mathematica script, and you literally have values you can compare to (not just "expect" based on physics/intuition), i guess that's a separate story. Maybe you can get more help if you present in pseudo-code form the high level problem, and ask specific matlab-related questions on how to accomplish a certain something that you are doing in Mathematica?