numerical solution of delay equation

matlabcode/delayedEquationModel/absDelayedMaturation.m

+classdef absDelayedMaturation < handle
+    properties (Access = protected)
+        timeI = [0:150]; %intergration time points
+        dataA; %data structure containing array of all data 
+        tpm;  %time of postmitotic start for non-core NC and Core region  [time_start_pm_NC time_start_pm_C]
+        ti; %time interphase assembly start  [time_start_ip_NC time_start_ip_C]
+    end
+    properties 
+        indir;
+        outdir;
+        dataF; %data structure containing data to be fitted
+    end
+    methods 
+         function MO = absDelayedMaturation(force, showplot)
+            if nargin < 1
+                force = 0;
+            end
+            if nargin < 2
+                showplot = 0;
+            end
+            MO.dataA = MO.dataIn(force, showplot);
+        end
+        function dataA = dataIn(MO, force, showplot)
+            %read in mature and pre
+            if nargin < 1
+                force = 0;
+            end
+            if nargin < 2
+                showplot = 0;
+            end
+            if ispc
+                indir = fullfile([getenv('HOMEDRIVE') getenv('HOMEPATH') ], 'Dropbox', 'NPCMaturation',...
+                    'expdata', 'emdata');
+                outdir = fullfile([getenv('HOMEDRIVE') getenv('HOMEPATH') ], 'Dropbox', 'NPCMaturation',...
+                    'results', 'emdata');
+            end
+            if ismac
+                indir = fullfile(getenv('HOME'), 'Dropbox', 'NPCMaturation',...
+                    'expdata', 'emdata');
+                outdir = fullfile(getenv('HOME'), 'Dropbox', 'NPCMaturation',...
+                    'results', 'emdata');
+            end
+            MO.indir = indir;
+            MO.outdir = outdir;
+            matfile = fullfile(indir,'emdata.mat');
+            if exist(matfile) && ~force
+                load(matfile)
+                return
+            end
+            
+            % first entry is minipore second entry is mature pore
+            data = struct('time', [], 'd_noncore', [], 'd_core_i', [], 'd_core_o', [], 't_core_i', [], 't_core_o', []);
+            dataA = [data; data];
+            xlsfile = fullfile(indir, 'Minipore-number-to-Antonio-141022.xlsx');
+            dataA(2).time = xlsread(xlsfile, 'Onlyvalues', 'C4:C15');
+            dataA(2).t_noncore = xlsread(xlsfile, 'Onlyvalues', 'D4:E15');
+            dataA(2).t_core_i = xlsread(xlsfile, 'Onlyvalues', 'F4:G15');
+            dataA(2).t_core_o = xlsread(xlsfile, 'Onlyvalues', 'H4:I15');
+            dataA(2).d_noncore = xlsread(xlsfile, 'Onlyvalues', 'D20:E31');
+            dataA(2).d_core_i = xlsread(xlsfile, 'Onlyvalues', 'F20:F31');
+            dataA(2).d_core_o = xlsread(xlsfile, 'Onlyvalues', 'H20:I31');
+            
+            dataA(1).time = xlsread(xlsfile, 'Onlyvalues', 'K4:K15');
+            dataA(1).t_noncore = xlsread(xlsfile, 'Onlyvalues', 'L4:M15');
+            dataA(1).t_core_i = xlsread(xlsfile, 'Onlyvalues', 'N4:O15');
+            dataA(1).t_core_o = xlsread(xlsfile, 'Onlyvalues', 'P4:Q15');
+            dataA(1).d_noncore = xlsread(xlsfile, 'Onlyvalues', 'L20:M31');
+            dataA(1).d_core_i = xlsread(xlsfile, 'Onlyvalues', 'N20:O31');
+            dataA(1).d_core_o = xlsread(xlsfile, 'Onlyvalues', 'P20:Q31');
+            
+            %% read and create surface area
+            surfdatafile = fullfile(indir, 'surfaceNucleus_24h_new.txt');
+            surfdata = load(surfdatafile);
+            dataA(1).surf = [surfdata(:,1)*60  surfdata(:,2:3)];
+            dataA(1).surf = [0 surfdata(1,2:3); dataA(1).surf];
+           
+            % perform a spline fit to smooth the data
+            dataA(1).surfspl = surfSpline( dataA(1).surf(:,1),  dataA(1).surf(:,2),7, 3, 8);
+            dataA(1).dsurfspl = fnder(dataA(1).surfspl, 1);
+            save(matfile, 'dataA');
+        end
+    end
+end

matlabcode/delayedEquationModel/ddex1.m

+function ddex1
+%DDEX1  Example 1 for DDE23.
+%   This is a simple example of Wille' and Baker that illustrates the
+%   straightforward formulation, computation, and plotting of the solution 
+%   of a system of delay differential equations (DDEs). 
+%
+%   The differential equations
+%
+%        y'_1(t) = y_1(t-1)  
+%        y'_2(t) = y_1(t-1)+y_2(t-0.2)
+%        y'_3(t) = y_2(t)
+%
+%   are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for
+%   t <= 0. 
+%
+%   The lags are specified as a vector [1, 0.2], the delay differential
+%   equations are coded in the subfunction DDEX1DE, and the history is
+%   evaluated by the function DDEX1HIST. Because the history is constant it
+%   could be supplied as a vector:
+%       sol = dde23(@ddex1de,[1, 0.2],ones(3,1),[0, 5]);
+%
+%   See also DDE23, FUNCTION_HANDLE.
+
+%   Jacek Kierzenka, Lawrence F. Shampine and Skip Thompson
+%   Copyright 1984-2004 The MathWorks, Inc.
+
+sol1 = dde23(@ddex1de,[40],@ddex1hist,[0 120]);
+sol2 = dde23(@ddex1de,[45],@ddex1hist,[0 120]);
+sol3 = dde23(@ddex1de,[50],@ddex1hist,[0 120]);
+sol4 = dde23(@ddex1de,[55],@ddex1hist,[0 120]);
+tint = [0:0.1:120];
+y = (deval(sol1,tint)+deval(sol2,tint)+deval(sol3,tint)+deval(sol4,tint))/4;
+figure;
+A=areaNuc(tint);
+plot(tint,y./repmat(A,3,1))
+title('An example of Wille'' and Baker.'); 
+xlabel('time t');
+ylabel('solution y');
+
+% --------------------------------------------------------------------------
+
+function s = ddex1hist(t)
+% Constant history function for DDEX1.
+    s = [0, 0, 2500];
+
+
+
+function [A, dAdt, kg] = areaNuc(t)
+par.a0 = 4.240e+02; 	 
+par.a1 = 1.613e+02; 	 
+par.kg1 = 7.219e-02 ;
+par.kg2 = 3.969e-01;
+    A = par.a0 + par.a1*(1-exp(-par.kg1*t)) + par.kg2*(t);
+    dAdt = par.a1*par.kg1*exp(-par.kg1*t)+par.kg2;
+    kg = (par.a1*par.kg1*exp(-par.kg1*t)+par.kg2)./A;
+
+
+% --------------------------------------------------------------------------
+
+function dydt = ddex1de(t,y,Z)
+% Differential equations function for DDEX1.
+ylag = Z(:,1);
+kM = 1.29;
+if t>10
+    v0 = 0.015;
+    v1 = 0.1;
+else
+    v0 = 0;
+    v1 = 0;
+end
+k = 0.2;
+kd = 4.2e-04;
+A = areaNuc(t);
+rho0 = 11;
+%dydt = [ (v0+v1*exp(-k*t))*A-kM*y(1)
+%         kM*y(1)-kM*ylag(1)
+%         kM*ylag(1)-kd*y(3)];
+dydt = [ v1*(A*rho0 - y(1)-y(2)-y(3))-kM*y(1)
+         kM*y(1)-kM*ylag(1)
+         kM*ylag(1)-kd*y(3)];

matlabcode/delayedEquationModel/delayedMaturationCls.m

+classdef delayedMaturationCls < absDelayedMaturation
+   methods 
+        function MO = delayedMaturationCls(force, showplot)
+            MO@absDelayedMaturation(force, showplot)
+            MO.dataF = MO.dataA;
+%             for i=1:2
+%                 MO.protF(3-i).core = MO.protA(i).tot_core;
+%                 MO.protF(3-i).noncore = MO.protA(i).tot_noncore;
+%             end
+        end
+        function so = solve(MO, par)
+            area = fnval(MO.dataF(1).surfspl,0);
+            
+            so = ode45(@MO.sumPores, [0:1:200], [area*4;area*6] , odeset, par, MO.dataF(1).surfspl)
+            
+        end
+        function dydt = sumPores(MO, t, y, par, surfspl)
+            area = fnval(surfspl,t);
+            
+            dydt(1) = par(1)*(par(2)*area - y(1)-y(2))-par(3)*y(1);
+            dydt(2) = par(3)*y(1);
+            dydt = dydt';
+            
+        end
+        
+        function plotResult(MO,par)
+            so = solve(MO, par);
+            plot(so.x, [so.y(1,:)./fnval(MO.dataF(1).surfspl, so.x); so.y(2,:)./fnval(MO.dataF(1).surfspl, so.x)])
+            %plot(so.x, par(1)*(par(2)*fnval(MO.dataF(1).surfspl, so.x)-so.y)./fnval(MO.dataF(1).surfspl, so.x))
+            %plot(so.x, so.y)
+        end
+            
+   end
+end
\ No newline at end of file

matlabcode/delayedEquationModel/delayedMaturation_miniPoreProduction.m

 %%
 % model of delayedMaturation with area growth using delayed equation
-%%
+% 
 function [parfit] = delayedMaturation_miniPoreProduction(par, timeShift, density)
 % interphase production Dultz and Ellenberg 2010 JCB
 % NRK cells 60-80 NPC/h however there is a dilution due to nuclear growth
@@ -56,7 +56,7 @@ par.ts = parArea(5);
 par.v1 = 0;     %rate of area dependent growth
 par.v2 = 0.014; %set to yield betweeh 180 min and 1200 min on average ~10.69 pores 
 par.v = par.v2;
-par.v3 = 0;
+par.v3 = 0; %a rate of post mitotic production. This is set to 0
 par.tShift = timeShift;
 par.P01 = 3;    %initial value Precursor in inner-core region 1 (at t=par.ts)
 par.P02 = 4;    %initial value Precursor in inner-core region2 (at t=par.ts)
@@ -394,7 +394,8 @@ end
 
 
 %%
-% compute distance to model
+% compute distance to model. Compute number of pores and then divide by
+% Area. This is much simpler then solving the equation for the denisties
 function [dist, dist2] = distModel(parfit, tofit_names, par, data, options)
 fitMeanMature = 0;
 %%
@@ -510,11 +511,11 @@ function NM = NMSol(t,y, par)
     %differential solution   
     [A, dAdt, kg] = areaNuc(t, par);
     if t<= par.tauM+par.tShift
-    NM = par.M0*exp(-par.kd*t);
+        NM = par.M0*exp(-par.kd*t);
     else
         C = 1/exp(-par.kM*t0)*(Np0 - par.v2*(par.a0/par.kM + par.a1*(1/par.kM - exp(-par.kg1*t0)/(par.kM - par.kg1)) + par.kg2*(t0/par.kM - 1/par.kM^2)) ...
         - par.v1*(par.kg1*par.a1*exp(-par.kg1*t0)/(par.kM-par.kg1)+par.kg2/par.kM));
-    NpInh = par.v2*(par.a0/par.kM + par.a1*(1/par.kM - exp(-par.kg1*t)/(par.kM-par.kg1)) + par.kg2*(t/par.kM-1/par.kM^2)) +...
+        NpInh = par.v2*(par.a0/par.kM + par.a1*(1/par.kM - exp(-par.kg1*t)/(par.kM-par.kg1)) + par.kg2*(t/par.kM-1/par.kM^2)) +...
         par.v1*(par.kg1*par.a1/(par.kM-par.kg1)*exp(-par.kg1*t) + par.kg2/par.kM);
         Np = NpInh + C*exp(-par.kM*t);
     end

matlabcode/odeModel/poreMaturationDen.m

@@ -91,7 +91,8 @@ classdef poreMaturationDen < absPoreMaturation
        function [idxm, idxd]  = findIdxs(MO)
            idxm = [1 1];
            idxd = [1 1];
-        end
+       end
+        
         function A = modelMatrix(MO, par)
         %% Generate matrix for multi step maturation process N1 steps with rate equal k1 and N2 steps with rate equal k2
         % k1 set the rate for formation of protein1 containing complex