- Author:
- Shelley Fong <sfon036@UoA.auckland.ac.nz>
- Date:
- 2022-06-10 15:16:24+12:00
- Desc:
- Move dSL to environment
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/7fb/rawfile/4e8b0e008b26cf300baf7a4c9f34903e8cfffaa8/Land_model_matlab/release/generate_figures.m
function generate_figures
global model
model = @Land2016;
addpath('./export_fig','./convert','./velocity_data');
set(0,'defaultlinelinewidth',1.5);
set(0,'DefaultAxesFontSize',10);
set(0,'DefaultTextFontSize',10);
humanCai = [0.166 0.166 0.165 0.165 0.164 0.164 0.165 0.166 0.170 0.175 0.181 0.188 0.196 0.205 0.215 0.225 0.235 0.247 0.259 0.270 0.280 0.289 0.297 0.304 0.312 0.319 0.328 0.336 0.345 0.353 0.359 0.364 0.368 0.373 0.377 0.382 0.388 0.394 0.401 0.408 0.413 0.416 0.420 0.423 0.427 0.432 0.438 0.444 0.450 0.455 0.460 0.463 0.465 0.468 0.470 0.475 0.480 0.485 0.491 0.496 0.499 0.502 0.504 0.504 0.506 0.510 0.513 0.519 0.524 0.529 0.531 0.532 0.532 0.532 0.533 0.535 0.538 0.543 0.547 0.551 0.554 0.554 0.554 0.554 0.553 0.555 0.559 0.561 0.565 0.568 0.570 0.569 0.568 0.568 0.566 0.567 0.570 0.574 0.577 0.581 0.581 0.581 0.579 0.578 0.578 0.578 0.579 0.582 0.586 0.588 0.588 0.587 0.585 0.583 0.581 0.580 0.583 0.587 0.590 0.593 0.593 0.592 0.590 0.586 0.584 0.585 0.586 0.590 0.592 0.594 0.593 0.593 0.590 0.588 0.586 0.586 0.587 0.589 0.592 0.594 0.593 0.591 0.589 0.586 0.584 0.583 0.585 0.587 0.589 0.591 0.589 0.588 0.585 0.582 0.580 0.579 0.580 0.582 0.585 0.587 0.586 0.584 0.581 0.579 0.578 0.577 0.578 0.580 0.581 0.583 0.581 0.579 0.576 0.572 0.570 0.569 0.572 0.574 0.576 0.578 0.577 0.574 0.572 0.569 0.567 0.566 0.567 0.569 0.570 0.570 0.569 0.567 0.563 0.560 0.558 0.557 0.558 0.559 0.561 0.562 0.563 0.560 0.557 0.554 0.551 0.550 0.550 0.552 0.553 0.555 0.553 0.550 0.546 0.544 0.541 0.541 0.542 0.543 0.545 0.545 0.543 0.541 0.537 0.534 0.531 0.531 0.530 0.532 0.534 0.535 0.533 0.530 0.527 0.523 0.520 0.518 0.518 0.519 0.521 0.521 0.520 0.518 0.515 0.512 0.509 0.508 0.508 0.508 0.510 0.510 0.508 0.505 0.502 0.499 0.497 0.496 0.496 0.496 0.496 0.496 0.496 0.493 0.490 0.486 0.484 0.483 0.483 0.483 0.484 0.483 0.481 0.479 0.476 0.472 0.469 0.467 0.467 0.467 0.468 0.469 0.466 0.464 0.460 0.457 0.454 0.452 0.453 0.453 0.454 0.454 0.453 0.450 0.447 0.444 0.441 0.439 0.439 0.440 0.441 0.441 0.440 0.437 0.433 0.430 0.427 0.426 0.426 0.426 0.426 0.426 0.424 0.421 0.418 0.415 0.413 0.411 0.410 0.411 0.412 0.411 0.409 0.407 0.404 0.401 0.399 0.398 0.397 0.397 0.397 0.397 0.396 0.393 0.390 0.387 0.384 0.383 0.383 0.383 0.384 0.384 0.382 0.380 0.377 0.374 0.372 0.370 0.369 0.370 0.370 0.370 0.368 0.366 0.364 0.361 0.358 0.357 0.357 0.357 0.357 0.356 0.355 0.353 0.350 0.348 0.345 0.344 0.343 0.344 0.344 0.344 0.343 0.341 0.338 0.336 0.333 0.332 0.332 0.332 0.332 0.332 0.332 0.330 0.327 0.325 0.323 0.322 0.321 0.321 0.321 0.321 0.320 0.318 0.316 0.313 0.311 0.310 0.310 0.310 0.311 0.311 0.310 0.308 0.306 0.303 0.301 0.301 0.300 0.300 0.301 0.301 0.300 0.298 0.296 0.294 0.292 0.291 0.291 0.291 0.292 0.292 0.291 0.289 0.287 0.285 0.284 0.283 0.282 0.283 0.283 0.282 0.281 0.280 0.279 0.277 0.275 0.274 0.274 0.274 0.275 0.274 0.273 0.272 0.270 0.269 0.268 0.267 0.267 0.267 0.267 0.267 0.266 0.265 0.263 0.261 0.260 0.260 0.260 0.260 0.260 0.260 0.259 0.258 0.257 0.255 0.253 0.253 0.253 0.253 0.254 0.254 0.253 0.252 0.250 0.248 0.247 0.246 0.246 0.247 0.247 0.247 0.247 0.246 0.244 0.243 0.241 0.241 0.241 0.242 0.242 0.242 0.242 0.240 0.239 0.237 0.236 0.236 0.236 0.236 0.237 0.237 0.236 0.236 0.234 0.233 0.232 0.231 0.232 0.232 0.233 0.232 0.232 0.231 0.229 0.228 0.227 0.226 0.226 0.227 0.228 0.228 0.227 0.226 0.225 0.224 0.223 0.222 0.222 0.223 0.223 0.223 0.223 0.222 0.221 0.220 0.219 0.218 0.219 0.219 0.220 0.220 0.220 0.219 0.218 0.217 0.216 0.215 0.215 0.216 0.216 0.216 0.216 0.215 0.214 0.213 0.212 0.211 0.212 0.212 0.213 0.213 0.213 0.212 0.211 0.209 0.209 0.209 0.209 0.209 0.210 0.210 0.209 0.208 0.207 0.206 0.205 0.205 0.205 0.206 0.206 0.206 0.206 0.206 0.205 0.204 0.203 0.203 0.203 0.204 0.204 0.204 0.204 0.203 0.202 0.201 0.200 0.200 0.200 0.201 0.201 0.201 0.201 0.200 0.199 0.199 0.198 0.197 0.198 0.198 0.199 0.199 0.199 0.198 0.197 0.196 0.195 0.195 0.196 0.196 0.197 0.197 0.197 0.196 0.195 0.194 0.194 0.193 0.194 0.194 0.195 0.195 0.195 0.194 0.193 0.192 0.191 0.191 0.191 0.192 0.192 0.193 0.192 0.192 0.191 0.190 0.190 0.189 0.190 0.190 0.191 0.191 0.191 0.190 0.189 0.189 0.188 0.188 0.188 0.188 0.189 0.189 0.189 0.188 0.187 0.187 0.186 0.186 0.186 0.187 0.187 0.188 0.187 0.187 0.186 0.185 0.185 0.184 0.185 0.185 0.186 0.186 0.186 0.185 0.184 0.184 0.183 0.183 0.183 0.184 0.185 0.185 0.185 0.184 0.183 0.182 0.182 0.181 0.182 0.182 0.183 0.184 0.184 0.183 0.182 0.181 0.181 0.181 0.181 0.181 0.182 0.182 0.182 0.182 0.181 0.180 0.179 0.179 0.180 0.180 0.181 0.181 0.181 0.180 0.180 0.179 0.178 0.178 0.178 0.179 0.180 0.180 0.180 0.179 0.179 0.178 0.177 0.177 0.177 0.178 0.179 0.179 0.179 0.178 0.178 0.177 0.176 0.176 0.176 0.177 0.178 0.178 0.178 0.177 0.176 0.175 0.175 0.175 0.175 0.176 0.176 0.177 0.177 0.177 0.176 0.175 0.175 0.175 0.175 0.176 0.176 0.177 0.176 0.176 0.175 0.174 0.174 0.174 0.174 0.174 0.175 0.175 0.175 0.175 0.174 0.173 0.173 0.173 0.173 0.174 0.174 0.174 0.174 0.174 0.173 0.172 0.172 0.172 0.172 0.172 0.173 0.173 0.173 0.173 0.172 0.171 0.171 0.171 0.171 0.172 0.172 0.173 0.173 0.172 0.171 0.171 0.170 0.170 0.170 0.171 0.172 0.172 0.172 0.172 0.171 0.170 0.170 0.169 0.170 0.171 0.171 0.171 0.171 0.171 0.170 0.169 0.169 0.169 0.169 0.170 0.170 0.171 0.171 0.170 0.170 0.169 0.168 0.168 0.168 0.169 0.170 0.170 0.170 0.170 0.169 0.168 0.168 0.167 0.168 0.169 0.169 0.170 0.170 0.169 0.168 0.168 0.167 0.167 0.168 0.168 0.169 0.169 0.169 0.169 0.168 0.167 0.167 0.167 0.167 0.168 0.168 0.169 0.169 0.168 0.167 0.167 0.166 0.167 0.167 0.168 0.168 0.168 0.168 0.168 0.167 0.166 0.166 0.166 0.166 0.166 0.167 0.168 0.167 0.167 0.166 0.166 0.165 0.165 0.166 0.167 0.167 0.168 0.168 0.167 0.166 0.166 0.165 0.165 0.166 0.166 0.167 0.167 0.167 0.167 0.166 0.165 0.165 0.165 0.166 0.166 0.167 0.167 0.167 0.166 0.166 0.165 0.165 0.165 0.165 0.166 0.166 0.167 0.167 0.166 0.166 0.165 0.165 0.165 0.165 0.166 0.166 0.167 0.167 0.166 0.165 0.165 0.164 0.164 0.165 0.166 0.166 0.166 0.166 0.166 0.165 0.165 0.164 0.164 0.165 0.165 0.166 0.167 0.167 0.166 0.165 0.165 0.164 0.164 0.165 0.165 0.166 0.166 0.166 0.166 0.165 0.165 0.164 0.164 0.165 0.166 0.166 0.166 0.166 0.166 0.165 0.165 0.164 0.164 0.164 0.165 0.166 0.166 0.166 ];
t_ca = 0:1000;
[t,y] = ode15s(@(t,y) model(t,y,humanCai,'intact',1,0), 0:1:(10*1000), model() );
tp('Model stats:',y(end-1000:end,1)*480);
lambdas = [0.9 0.95 1 1.05 1.1];
figure(101);
clf;
subplot(1,2,2); hold on;
for l = lambdas
[t,y] = ode15s(@(t,y) model(t,y,humanCai,'intact',l,0), 0:1:(10*1000), model() );
if l==1; lw=2; else lw=1; end
plot(t-9e3,y(:,1)*480,'LineWidth',lw);
end
xlabel('Time (ms)');
ylabel('Force (kPa)');
legend( arrayfun(@(l) sprintf('\\lambda = %g',l),lambdas,'UniformOutput',0));
xlim([0 1000]);
set(gca,'XTick',[0 250 500 750 1000],'YTick',[0 25 50 75 100]);
subplot(1,2,1); hold on;
plot(t_ca,humanCai);
xlabel('Time (ms)');
ylabel('Intracellular calcium concentration (\mu{}M)');
ylim([0 0.6]);
set(gca,'XTick',[0 250 500 750 1000],'YTick',[0 0.15 0.3 0.45 0.6]);
myfig2(figure(101), './Fig_fit_ca_twitch.png',[7 2.5]*2,225);
% --- velocity
figure(102); clf;
t={}; tension={}; traces={};
[t{1},tension{1},traces{1}] = plot_data(0);
[t{2},tension{2},traces{2}] = plot_data(1);
[t{3},tension{3},traces{3}] = plot_data(2);
Cai = 30;
[~,y] = ode15s(@(t,y) model(t,y,Cai,'skinned',1,0),0:10000,model());
y00 = y;
[score, tension_out,rel] = score_par(y00,t,tension);
figure(102); clf;
for i=1:3
subplot(1,4,i); hold on;
plot(t{i},tension{i}/30,'Color',[ 0 0.4470 0.7410],'LineWidth',2);
plot(t{i}(1:end-1),tension_out{i}/30,'Color',[ 0.8500 0.3250 0.0980],'LineWidth',2);
for j=1:length(traces{i})
plot(t{i},traces{i}{j}/30,'Color',[0.05 0.05 0.05]*(j+4),'LineWidth',1);
end
plot(t{i},tension{i}/30,'Color',[ 0 0.4470 0.7410],'LineWidth',1);
plot(t{i}(1:end-1),tension_out{i}/30,'Color',[ 0.8500 0.3250 0.0980],'LineWidth',2);
end
legend('Mean data','Model fit','Individual traces');
subplot(1,4,1); ylabel('Normalized force');
for i=1:3
subplot(1,4,i);
xlabel('Time (ms)');
xlim([-75 2000]); set(gca,'XTick',[0 500 1000 1500 2000],'YTick',[0.5 0.75 1 1.25 1.5]);
ylim([0.5 1.6]);
end
%myfig2(figure(102), './Fig_fit_vel.png',[7 2.5]*2,225);
% --- constant shortening
files = {'constant_shortening/16318001_corrd.abf','constant_shortening/16318002_corrd.abf','constant_shortening/16318003_corrd.abf','constant_shortening/16318004_corrd.abf','constant_shortening/16318005_corrd.abf','constant_shortening/16318006_corrd.abf','constant_shortening/16318007_corrd.abf','constant_shortening/16318008_corrd.abf','constant_shortening/16318009_corrd.abf','constant_shortening/16318010_corrd.abf','constant_shortening/16318011_corrd.abf'};
area = [52*13 52*13 52*13 52*13 52*13 52*13 29*14 29*14 29*14 29*14 29*14];
sl0 = [2.1 2.1 2.06 2.1 2.3*ones(1,11)];
slrest = 1.9;
length_change = [8.5 8.5 8.5 10.0 11.4 11.4 10.0 8.5 7.1 7.1 8.5]' * [1 1 1 1]; % regular, 1-11
length_dur = [ 100 200 300 0
100 200 300 0
150 200 250 300
150 200 250 300
150 200 250 300
100 200 300 0
100 200 300 0
100 200 300 0
100 200 300 0 % 009
150 200 250 300
150 200 250 300
150 200 250 300 % 011
];
t0 = 8310;
tsh = 11000;
trx_sh = 14000;
trx_long = 19500;
x = 0;
lambdas = {};
fnorms = {};
ts = {}; tend={};
lx=[]; fe=[];
for i=1:length(files)
[d,si,h] = abf2load(files{i});
t = (0:size(d,1)-1)*(si/1000);
f = 1000 * squeeze(d(:,1,:)) / area(i);
l = squeeze(d(:,2,:));
x0 = x+1;
fni = [];
for j=1:size(f,2)
x=x+1;
ts{x} = t;
[~,ti] = min(abs(t-(t0+length_dur(i,j))));
[~,tsi] = min(abs(t-t0));
tend{x} = t(ti);
fprintf('%d/%d: mean f = %f\n',i,j,mean(f(t < t0-5 & t > t0-40,j)));
frx_sh(x) = mean( f( abs(t-trx_sh)<50 ,j) );
frx_long(x) = mean( f( abs(t-trx_long)<50 ,j) );
fsh(x) = mean( f( abs(t-tsh)<50 ,j) );
fe(x) = mean( f(ti:ti+100,j)); % ~ 20ms
fs(x) = mean( f(tsi-100:tsi,j)); % ~ 20ms
cor_fe(x) = (fe(x) - frx_sh(x)) ./ (fs(x) - frx_long(x));
dur(x) = length_dur(i,j); lx(x) = length_change(i,j)/100;
end
end
figure(102); subplot(1,4,4)
plot(lx./dur,cor_fe,'o', 0.1./(100:10:300), rel(:,1)./rel(:,2), 'x-',[0.2/200 0.2/150 0.2/150 0.2/200 0.2/200],[0.9 0.9 0.5 0.5 0.9],'--'); % [100 200 300]
ylim([0 0.91]);
xlabel('Shortening velocity (/ms)'); ylabel('Relative force at end of shortening'); legend('Experimental data','Skinned Model prediction','Expected force in organ','Location','NorthWest')
ax1 = gca;
ax1_pos = ax1.Position; % position of first axes
ax2 = axes('Position',ax1_pos,...
'XAxisLocation','top',...
'YAxisLocation','left',...
'Color','none');
xlabel('Sarcomere length shortening speed (\mu{}m/s)')
scale = 2.3*1000;
ax2.XLim = ax1.XLim*scale;
ax2.XTick=ax1.XTick*scale
ax2.YTick=[];
myfig2(figure(102), './Fig_fit_vel_cmb.png',[10 2.5]*2,225);
function [tpt,rt50,rt95,mif,mxf] = tp(tag,fi)
[mxf,tpt] = max(fi);
f0 = fi-min(fi);
f0(1:tpt) = f0(tpt);
t50 = find(f0 <= 0.5 * max(f0),1);
t90 = find(f0 <= 0.1 * max(f0),1);
t95 = find(f0 <= 0.05 * max(f0),1);
if isempty(t95); t95=1000; end
onset = find(diff(fi)>0.01,1); % ?
fprintf('%s\tTPT: %.1f\tRT50: %.1f\tRT90: %.1f\tRT95: %.1f\tMax %.1f\tMin %.3f\n',tag,tpt-onset,t50-tpt,t90-tpt,t95-tpt,mxf,min(fi));
rt50 = t50-tpt;
rt95 = t95-tpt;
tpt = tpt-onset;
mif = min(fi);
function [score,tension_out,rel,allrel] = score_par(y00,t,tension)
global model
stretch = [0.5 1 2];
step_dur = 10;
step_t = 2000;
Cai = 30;
xb0 = [0.95993 0.97022 0.9748];
ldp = [ 5.25 6.6 4.5;
-2.4 -2.4 -2.4];
tscale= [0.8498 0.8396 0.8388];
dif = [0 0 0];
for i=1:3
Tfit = tension{i};
flen = ones(size(Tfit));
flen(step_dur:step_t) = stretch(i)/100 + 1; % stretched
st = interp1([0 step_dur],[1 flen(6+step_dur)],0:step_dur);
flen(1:step_dur+1) = st;
flen(step_t:step_t+step_dur) = st(end:-1:1);
ti = t{i};
f = @(tt,y) model(tt,y,Cai,'skinned',interp1(ti,flen,tt), (interp1(ti,flen,tt+1e-3)-interp1(ti,flen,tt))/1e-3 , ldp(:,i)); % TODO: LD
y0 = y00( find(y00(:,1)/y00(end,1) < xb0(i), 1,'last'), :) %
[tt,y] = ode15s(f,t{i}(1:end-1),y0,odeset('MaxStep',1)); % reduce xb0
tatp = cell2mat(cellfun(@(ty) get_ta(f,ty),num2cell([tt y],2),'UniformOutput',0));
ta = tatp(:,1);
tp = tatp(:,2);
tot_tension = (ta + tp)*tscale(i);
tension_out{i} = tot_tension;
[mt(i),mi] = max(tension{i});
mm(i) = max(tot_tension);
mm0(i) = tot_tension(1);
dif(i) = norm(tot_tension - tension{i}(1:end-1));
end
q=1; i=1;
goal = [0.2 0.4 0.5];
for te = 100:10:300
f2 = @(tt,y) model(tt,y,Cai,'skinned', flen(1) * (1 - (0.1/te)*tt) , -0.1/te , ldp(:,end));
[tt2,y2] = ode15s(f2,0:te,y0,odeset('MaxStep',1)); % reduce xb0
tatp2 = cell2mat(cellfun(@(ty) get_ta(f,ty),num2cell([tt2 y2],2),'UniformOutput',0));
ta2 = tatp2(:,1);
rel(i,:) = [ta2(end) ta2(1)]; i=i+1;
if te==100||te==200||te==300
shs(q) = 50*abs( (ta2(end)/ta2(1)) - goal(q) );
q=q+1;
end
end
fprintf('score = %f + %f\n',sum(dif),sum(abs(shs)))
score = sum(dif) + sum(abs(shs));
function [t,ta,tp]=get_ta(model,ty)
[~,~,ta,tp] = model(ty(1),ty(2:end));
t = [ta tp];
