How to fit data to a PDE using lsqcurvefit/lsqnonlin where y_pdepe is calculated from pdepe output

Hi All,

I want to fit some experimental data to my system of PDEs and estimate some system paramters. Trouble is my numerical output is not the direct output of my PDE rather calculated from the PDEPE output (i.e. I am caculating the flux from the PDEPE calculated conc. profile and then using that to calculate current). This is because of my experimental data (current flux vs substrate concentration).

}

Now my experiemntal data is in the form

s_bulk = [0, 5e-06, 10e-06, 20e-06, 50e-06, 100e-06]

I_exp = [0, 2.35e-03, 2.7e-03, 3.125e-03, 3.29e-03, 3.25e-03]

What I want to do is to fit my pdepe output "I_num" to the "I_exp" given above.

My lsqcurvefit call:

% The purpose of this script is to call PSw_EK_v8 using multiple
% substrate concentrations and calculate the respective current values.
% Using these values we can perhaps try to find the k_m value.  

% k_m            = [4.92e07, 4.92e07];
s_bulk         = [0, 5e-06];                                          % mol/cm^3 
I_exp          = [0, 2.35e-03];                                       % mA/cm^2
% I_num          = arrayfun(@PSw_EK_v8, k_m, s_bulk, 'un', 0);
% I_num          = cell2mat(I_num);

k_m_0          = 4.92e06;
k_m            = lsqcurvefit(@PSw_EK_v8, k_m_0, s_bulk, I_exp);

My pdepe function is given as:

function [I_num] = PSw_EK_v8(k_m, s_bulk)

C.n         = 1;                                 % No. of Electrons Involved
C.F         = 96485.3329;                        % C/mol
C.R         = 8.314;                             % J/(mol.K) or CV/(mol.K)
C.Temp      = 273.15+37;                         % K

C.D_s       = 7.00E-06;                          % cm^2/s
C.D_m       = 2.81E-06;                          % cm^2/s (De of polymer)
C.D_e       = 0;                                 % because covalently bound

% C.k_m       = 4.92E+07;                          % (mol/cm3)^{-1}s^{-1}
C.kcat      = 800;                               % Turnover No. (1/s) 
C.k_e       = C.kcat;                            % Refer to Schnell
C.KM        = 2.00E-05;                          % Michalis Constant (mol/cm^3)
C.P_s       = 1;                                 % Partition Coefficient
C.P_m       = 1;                                 % Partition Coefficient

% Initial Concentrations 
C.e_layer   = 5.00E-08;                          % mol/cm^3
C.m_layer   = 1.3E-06;                           % mol/cm^3

C.Area      = 0.0707;                            % cm^2 
C.l         = 0.001;                             % cm 
m           = 0;                                 % Cartesian Co-ordinates

N           = 30;                                % No. of Points
 
C.tspan1    = linspace(0, 500, N);               % s 
% C.tspan1    = logspace(log10(0.000006), log10(500), N); C.tspan1(1) = 0;

xmesh       = linspace(0, C.l, N);               % cm 
% xmesh       = logspace(log10(0.000006), log10(C.l), N); xmesh(1) = 0;

% FIRST LINEAR POTENTIAL SWEEP PDEPE Solver Command
E           = 0.0;                               % V
C.E0        = 0.2;                               % V
C.ScanRt    = 0.001;                             % V/s
C.E_1       = E+(C.ScanRt*C.tspan1);             % Potential Sweep 1
C.epsilon1  = ((C.n*C.F)/(C.R*C.Temp)).*(C.E_1-C.E0); 
C.c_mred1   = C.m_layer./(1+exp(C.epsilon1));    % Mred

% FIRST POTENTIAL SWEEP Initial & Boundary Condition Call
ic_arg      = {@(x)s_bulk*ones(size(N)) ; @(x)C.e_layer*ones(size(N)); ...
               @(x)C.m_layer*ones(size(N))}; 
IC          =  @(x)PDE_PSw_EK_IC(x, ic_arg); 
BC          =  @(xl, cl, xr, cr, t)PDE_PSw_EK_BC(xl, cl, xr, cr, t, C, s_bulk); 

% Calling pdepe solver after time mesh refinement 
% MaxStep sets an upper bound on the size of any step taken by the solver. 
% If the equation has periodic behavior, for example, then setting MaxStep 
% to a fraction of the period ensures that the solver does not enlarge the 
% step so much that it steps over an area of interest.
optns      = odeset('MaxStep',1e-01,'RelTol',1e-09,'AbsTol',1e-10); 
sol1       = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, C, k_m), ... 
                   IC, BC, xmesh, C.tspan1, optns);
               
% Concentration Profiles c(i, j, k)(Solutions) 
c1          = sol1(:, :, 1);                     % Substrate Conc.  
c2          = sol1(:, :, 2);                     % Ered Conc.
c3          = sol1(:, :, 3);                     % Mred Conc.
c4          = C.m_layer-c3;                      % Mox Conc. 
   
% Plot Species Conc. in the Layer for FIRST POTENTIAL SWEEP 
figure(1); 
[hAx,hLine1,hLine2]= plotyy(xmesh,c4(N,:),xmesh,c2(N,:)); 
xlabel('x / cm'); 
ylabel(hAx(1),'[c_{m_{ox}}] / mol cm^{-^3}')     % left y-axis
ylabel(hAx(2),'[c_{e_{red}}] / mol cm^{-^3}')    % right y-axis 
title('Concentratrions for FIRST POTENTIAL SWEEP');

% % Plot Species Conc. in the Layer 
% figure(11)
% plot(xmesh, c1(N,:));
% xlabel('x / cm'); ylabel('substrate [s]');

%     [~, dc1Sdx_0(counter)]  = pdeval(m, xmesh, c1(counter,:),  0);  % Substrate Flux at x=0
%     [~, dc1Sdx_l(counter)]  = pdeval(m, xmesh, c1(counter,:),C.l);  % Substrate Flux at x=d
[~, dc4Mdx_0]  = pdeval(m, xmesh, c4(N,:),  0);  % Mox Flux at x=0
%     [~, dc4Mdx_l(counter)]  = pdeval(m, xmesh, c4(counter,:),C.l);  % Mox Flux at x=d

j_num_ano  = (-C.D_m .* dc4Mdx_0);
I_num      = (C.n * C.F .* j_num_ano)/(C.Area); 
% I vs E for FIRST POTENTIAL SWEEP 
figure(4); 
plot(C.E_1, I_num, 'r--', 'LineWidth',0.5);
xlabel('E / V'); ylabel('I / A'); 
hold on;  

z = 1; 

end

%% pdepe Function PSw_EK_v8
%%

% Notes C.D_e is 0 because enzyme is immobilized in the polymeric matrix 
% and cannot diffuse across the layer.  

function [cc, ff, ss] = PDE_PSw_EK(x, t, c, DcDx, C, k_m)

%          Substrate;  Ered;  Mred;   
cc          = [    1;     1;     1]; 
ff          = [C.D_s; C.D_e; C.D_m].*DcDx;

menten      = (C.k_e*c(1))*(C.e_layer-c(2))/(C.KM+c(1));
mred        = k_m*(C.m_layer-c(3))*c(2);

ss          = [-menten; -mred+menten; +mred]; 
           
end

%% Initial Condition --> Initial Concentrations of Species 
%%

function c0 = PDE_PSw_EK_IC(x, ic_arg)
   
%                 Substrate;         Ered;          Mred;   
c0          = [ic_arg{1}(x); ic_arg{2}(x); ic_arg{3}(x)];

end

%% Boundary Condition - 1 
%%

function [pl, ql, pr, qr] = PDE_PSw_EK_BC(xl, cl, xr, cr, t, C, s_bulk)

% ---------------------
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% -------|-------------
%        pl           pr 

% Order --> Substrate; Mox; Mred; Eox; Ered 

%              Substrate;Er;                                  Mox;   
pl     = [0             ; 0; cl(3)-interp1(C.tspan1,C.c_mred1,t)];  
 
ql     = [1             ; 1; 0];

pr     = [cr(1)-s_bulk  ; 0; 0];

qr     = [0             ; 1; 1]; 

end

[Star Strider]

I never tried to do this with pdepe (and I have very little experienced with pdepe) however a typical approach for ordinary differential equations is presented in Parameter Estimation for a System of Differential Equations and similar threads.

If the partial differential equations are linear and have constant coefficients (they appear to, however I could be missing something), another option is to use the Symbolic Math Toolbox to first take the Laplace transform of the time derivatives, creating a system of ordinary differential equations in the spatial (or other) derivatives that can then be solved analytically. Then invert them back to the time domain and fit the data with that. (That was a grad school final exam problem, and something of a challenge to do my hand, however I got it right.)

.

[Hashim]

I guess what I want to do here is to get multiple values of i_num against s_0 and then use lsqcurvefit to estimate k_m. Anyway you think that might be possible please let me know. I have been able to get multiple values of i_num by passing s_0 using arrayfun but do not know how to integrate that with lsqcurvefit so that it automatically gives me an estimate for k_m.

k_m = arrayfun(@PSw_EK_v8, s_bulk, 'un', 0)

Also, I have gone through pretty much all the threads you've mentioned.

For the second part of your suggestion do you mean to say that the equations provided can be solved analytically using symbolic toolbox? because my assumption so far was that the symbolic toolbox is likely to give an significantly long solution (as per an experienced mathematica user). If yes, could you point me towards some resources?