Hi, I am trying to implement the Electroporation model described in the paper "Model of Creation and Evolution of Stable Electropores for DNA Delivery"
It is basically a system of differential equations with discontinuities. I used an events function to handle the discontinuites, but it seems like either the events function is not working, or it is taking too long to run because of the large number of differential equations, or both.
Any suggestions?
function Modelclc; close all; clear all;D0 = 2*10^-12 ; % DNA diffusion coefficient
R = 10*10^-6; % cell radius
H = 5*10^-9; % membrane thickness
Vcell = 3.14*10^-16; % cell volume
C0 = 2*10^-6; % Outside DNA concentraion
zeff = -27; % effective valence of DNA
q0 = 1.602*10^-19; % elementary charge
kB = 1.381*10^-23; % Boltzman constant
T0 = 310; % Absolute Temperature (37 C)
g = 2; % Conductivity of the solution
f = 1.5; % form factor
dVs = 1; % Potential Threshold
R0 = 100; % external resistance (series resistance of experimental setup)
Rm = 0.523; % Surface resistance of the membrane
Cm = 9.5*10^-3; % Surface Capacitance of the membrane
F_max = 0.7*10^-9; % Maximum electric force for Vm = 1V
rh = 0.97*10^-9; % Constant
rt = 0.31*10^-9; % Constanta = 1*10^9; % creation rate coefficient
Vep = 0.258; % characteristic voltage of electroporation
N0 = 1.5*10^9; % equillibrium pore denisty at Vm = 0
rp = 0.51*10^-9 ; %minimum radius of hydrophillic pores
rm = 0.8*10^-9; % minimum energy radius at Vm = 0
D = 5*10^-14 ; % Difussion coefficient for pore radius
Beta = 1.4*10^-19; % Steric repulsion energy
gamma = 1.8*10^-11; %Edge Energy
s0 = 1*10^-3; %Tension of the bilayer without pores
s1 = 2*10^-2; %Tension of hydrocarbon-water interface
tau = 2.46; % electroporation constant
kf = 2*10^-9; %coefficient
k0 = 1*10^10; %coefficientq = (rm/rp)^2; % model paramater
A = 4*pi*R^2; % Total Surface Area
C = Cm*A; % total capacitance
R = Rm/A; % total resistance
Ip = 0 ; % initialize combined current through all pores
K = 0; % initialize number of large pores in a given itteration
K1 = 0; % initialize total number of large pores
% s_eff = s0; % initial s_eff
% A_p = 0; % initial A_p
% Runge-Kutta
% ------------------------------------------------------------------------
timerange = [0 1*10^-7]; % microseconds
IC = [N0;0;zeros(1000,1)]; % Initial Conditions : N(0) = N0, Vm(0) = 0
y_store = transpose(IC(1:2)); % storage for N and Vm
t_store = [ 0 ]; % storage for time
T = (10^-9)*(0:100) ; % Cooresponding Time vector for Applied Voltage
V0 = ones(1,101); % Applied Voltage 0.1 microsecond (100 microseconds total)
Pores = zeros(1000,1000); % Initialize Pore history storage array.
%ith row corresponds to ith itteration, jth collumn corresponds to jth pore.
options = odeset('events', @Efcn); % set options (events function)
for i = 1:1000 [t , y] = ode23(@(t,y) odes(t, y, T , V0) , timerange , IC, options); IC = y(size(y,1),:); % initial conditions for next timestep (last row of y)
IC = IC(:); % change to collumn vector
y_store = [y_store ; y(:,1:2)]; % stores values for N and Vm (m x 2)
t_store = [t_store ; t + t_store(length(t_store))]; % stores values (m x 1)
if K1 >= 1for k = 1:K1 % itterate through pores
Pores(i,k) = y(size(y,1),2+k); % store radius values
if y(size(y,1),2+k) < rm && y(size(y,1),2+k) ~= 0 % if pore radius has fallen below rm
Pores(i,k) = 0; % set pore radius to 0 indicating that it is absorbed into small pore population
IC(2+k) = 0; % set initial condition for corresponding pore to r = 0;
IC(1) = IC(1) + 1/A; % add one pore to small pore population
end end end K = 0; K = K + floor(y(size(y,1),1)*A); % Launch Pores
if K >= 1 IC(1) = IC(1) - floor(y(size(y,1),1)*A)/A; % subtract integer number of pores from small pore population
rng(0,'twister'); sd = 0.001*10^-9; % standard deviation 0.001 nm
mean = rm; % mean of rm
new_pores = sd.*randn(1,K) + mean; % create normal distribution of pores with mean = rm and sd = 0.001 nm
Pores(i+1,K1+1:K1+K) = new_pores; % save pores
end if K>=1IC(3+K1:2+K1+K) = new_pores(:) ; % add new pores to initial conditions
K1 = K1 + K; % Total pores
end assignin('base','pores',Pores)assignin('base','Vm', y_store(:,2))assignin('base','time', t_store)assignin('base','pore_density',y_store(:,1))end function rk = odes(t, y, T, V0)format longV0 = interp1(T, V0, t);rk(1) = a*exp((y(2)/Vep)^2)*(1 - (y(1)/(N0*exp(q*(y(2)/Vep)^2)))); % dN/dt
% Ip = (((y(2)/((H/(pi*g*(rm^2))) + (1/(2*g*(rm)))))*y(1)*A) + 2*pi*g*y(2)*sum((y(3:end).^2)./(2*pi + pi*y(3:end)),'all'));
rk(2) = ((V0/R0) - (((y(2)/((H/(pi*g*(rm^2))) + (1/(2*g*(rm)))))*y(1)*A) + 2*pi*g*y(2)*sum((y(3:end).^2)./(2*pi + pi*y(3:end)),'all')) - (((1/R0) + (1/R))*y(2)))/C; % dVm/dt
for j = 3:K1+2 %drj/dt
rk(j) = (-D/kB*T0)*((-4*Beta*rp*(y(j))^-5) + 2*pi*gamma - 2*pi*2*s1 - ((2*s1 - s0)/(1 - (((y(1))*A*pi*(rm^2) + (sum(((y(3:K1)).^2),'all'))*pi)/A)))*((y(1))*A*pi*(rm^2) + (sum(((y(3:K1)).^2),'all'))*pi)*y(j) + (F_max*((y(2))^2)/(1 + (rh/(y(j) + rt)))));end for v = 3:1002 if y(v) == 0 rk(v) = 0; % drj/dt = 0 for pores that have not been created, or that have been absorbed
end end rk = rk(:) ; end end function [check,isterminal,direction] = Efcn(t,y)isterminal = 1;direction = [];check = double(floor(y(1)*(4*pi*(10*10^-6)^2)) < 1); % checks if N*A is greater than 1
for m = 3:1002 if y(m) < 0.8*10^-9 && y(m) ~= 0 % Absorb pores that have fallen below rm
check = 0; end endend
Best Answer