close all clear clc % Code for creating the figures and processing the .mat variable containing % all the subject data % Written by L. Welte 2020 for : The elasticity of the plantar fascia influences the windlass effect on arch deformation during human running % Proceedings B % All methods are presented in the manuscript. % Note: all base figures were created in MATLAB then edited and combined % with other images in PowerPoint. subj_order = [12,10,2,4,1,3,5,6,7,8,9,11]; % indices 12, 10 and 2 are RFS, the rest are FFS, 12 has extra ligaments % IMPORTANT: subject indices do NOT correspond to the participant numbers % in the manuscript. They are decreasing from 12:1 in the order presented above. % Participant 12 in the manuscript is variable index 12, % participant 11 is index 10, participant 10 is index 2 etc. nsubj = 12; % number of subjects load('results.mat') % these must be on your MATLAB path load('colorblind_map.mat') %% Make figure 2 % The frames over which an average is stable frs_av = 1:90; frs_av_toe = 1:95; % 2a Arch dorsiflexion figure; hold on; patch([0 0 20 20 ], [30 -27 -27 30],'k','Facealpha',0.1,'Linestyle','none'); patch([55 55 85 85 ], [30 -27 -27 30],'k','Facealpha',0.1,'Linestyle','none'); h1 = plot( results_norm.angles_filt.cal_mt1.DP(subj_order(1:3),:)','k--'); h2 = plot( results_norm.angles_filt.cal_mt1.DP(subj_order(4:end),:)','k'); plot(frs_av-1,nanmean( results_norm.angles_filt.cal_mt1.DP(:,frs_av)),'k','Linewidth',4) xlim([0 100]) ylim([-27 30]) xlabel('% Stance') ylabel(sprintf('Arch Dorsiflexion [%s]',0176)) legend([h1(1),h2(1)],{'RFS','FFS'}) % 2b arch adduction figure; hold on; h1 = plot( results_norm.angles_filt.cal_mt1.AD(subj_order(1:3),:)','k--'); hold on; h2 = plot( results_norm.angles_filt.cal_mt1.AD(subj_order(4:end),:)','k'); plot(frs_av-1,nanmean( results_norm.angles_filt.cal_mt1.AD(:,frs_av)),'k','Linewidth',4) legend([h1(1),h2(1)],{'RFS','FFS'}) xlim([0 100]); xlabel('% Stance') ylabel(sprintf('Arch Adduction [%s]',0176)) %% Make figure 3 % 3a MTPJ dorsiflexion figure; hold on; patch([0 0 20 20 ], [-10 80 80 -10],'k','Facealpha',0.1,'Linestyle','none'); patch([55 55 85 85 ], [-10 80 80 -10],'k','Facealpha',0.1,'Linestyle','none'); h1 = plot( results_norm.angles_filt.mt1_ph1.DP(subj_order(1:3),:)','k--'); hold on; h2 = plot( results_norm.angles_filt.mt1_ph1.DP(subj_order(4:end),:)','k'); plot(frs_av_toe-1,nanmean( results_norm.angles_filt.mt1_ph1.DP(:,frs_av_toe)),'k','Linewidth',4) legend([h1(1),h2(1)],{'RFS','FFS'}) xlim([0 100]); xlabel('% Stance') ylabel(sprintf('MTPJ Dorsiflexion [%s]',0176)) % 3b plantar fascia elongation figure; hold on; patch([0 0 20 20 ], [30 -20 -20 30],'k','Facealpha',0.1,'Linestyle','none'); patch([55 55 85 85 ], [30 -20 -20 30],'k','Facealpha',0.1,'Linestyle','none'); h1 = plot( results_norm.lig_filt.pf_strain(subj_order(1:3),:)','k--'); hold on; h2 = plot(results_norm.lig_filt.pf_strain(subj_order(4:end),:)','k'); plot(frs_av-1,nanmean( results_norm.lig_filt.pf_strain(:,frs_av)),'k','Linewidth',4) legend([h1(1),h2(1)],{'RFS','FFS'}) xlim([0 100]); ylim([0.9,1]) xlabel('% Stance') ylabel(sprintf('Plantar Fascia Elongation ')) % 3c arch vel figure; hold on; patch([0 0 20 20 ], [0 900 900 0],'k','Facealpha',0.1,'Linestyle','none'); patch([55 55 85 85 ], [0 900 900 0],'k','Facealpha',0.1,'Linestyle','none'); h1 = plot( results_norm.angles_filt.arch_vel(subj_order(1:3),:)','k--'); hold on; h2 = plot(results_norm.angles_filt.arch_vel(subj_order(4:end),:)','k'); plot(frs_av-1,nanmean( results_norm.angles_filt.arch_vel(:,frs_av)),'k','Linewidth',4) legend([h1(1),h2(1)],{'RFS','FFS'}) xlim([0 100]); xlabel('% Stance') ylabel(sprintf('Magnitude of arch angular velocity [%s]',0176)) %% make figure 4 a thru c thrA = 0.5; % MTP threshold thrB = 0.0005; % PF threshold % make 4a-c arch_vel = mean(results_norm.angles_filt.arch_vel((subj_order),:)); varA = mean(results_norm.angles_filt.mt1_ph1.DP((subj_order),:)); varB = mean(results_norm.lig_filt.pf_strain((subj_order),:)); codingPFphases(varA,varB,thrA,thrB,arch_vel); % make 4d arch_vel = arch_vel+1000; % this just removes it from the third panel so that you can see all the participants varA = results_norm.angles_filt.mt1_ph1.DP((subj_order),:); varB = results_norm.lig_filt.pf_strain((subj_order),:); [codes_all,codes_toe,~] = codingPFphases(varA,varB,thrA,thrB,arch_vel); %% make Figure 5 s_count = 1; for subj = subj_order % get simple variables mtp_angle = results_norm.angles_filt.mt1_ph1.DP(subj,:); arch_def = results_norm.kinemat_filt.arch_def(subj,:); pf_length = results_norm.lig_filt.pf_strain(subj,:); % calculate the differences arch_diff = calculateVelocity(arch_def,1); mtp_diff = calculateVelocity(mtp_angle,1); % get rid of the values at the beginning that are not part of % propulsion ind_sm = false(1,101); ind_sm(1:40) = 1; mtp_diff(ind_sm) = NaN; arch_diff(ind_sm) = NaN; delta_s(s_count,:) = results_norm.r_sens(subj) * mtp_diff; % met head radius * pi/180 * delta MTP ; note r_sens has pi/180 built in mtp_diff_s(s_count,:) = mtp_diff; arch_diff_s(s_count,:) = arch_diff; s_count = s_count+1; end diff_btw_models = -1*(arch_diff_s+delta_s); % negative to match convention ind_inh_wind = codes_all == 1; ind_pure_wind = codes_all == 2; ind_enh_wind = codes_all == 3; ind_inh_wind(:,1:40) = 0; ind_pure_wind(:,1:40) = 0; ind_enh_wind(:,1:40) = 0; % initialize sens_inh = diff_btw_models ; sens_wind = diff_btw_models; sens_enh = diff_btw_models; % set all non points in each phase to NaNs sens_inh(~ind_inh_wind) = NaN; inh = nanmean(sens_inh,2); sens_wind(~ind_pure_wind) = NaN; wind = nanmean(sens_wind,2); sens_enh(~ind_enh_wind) = NaN; enh = nanmean(sens_enh,2); makeStatsDotPlot(1:3,[inh,wind,enh]',{'inh','wind','enh'},cmap_cb([8 3 5],:),'diamond') plot_pts = [inh,wind,enh]'; for subj = 1:nsubj hp = plot(1:3,plot_pts(:,subj),'marker','none','color','k'); if subj < 4 set(hp,'Linestyle','--') end end ylabel('Arch length change relative to MTPJ arc length [mm]') %% Make figure 6 ligs = {'PLPL','PCCL','PCNL_inf','PCNL_med','pf'}; % list of ligaments norm2minmax = @(x) (x-min(x))/(max(x)-min(x)); % function to normalize between 0 and 1 fig = figure; set(fig,'defaultAxesColorOrder',[ [159, 131, 199]/256;[0 0 0]]); % set up the elongation yyaxis left hold on; for L = 1:length(ligs) h(L)=plot(norm2minmax(results_norm.lig_filt.(ligs{L})(12,:))); end h(5).Color = [79, 177, 196]/256; ylabel('Normalized Elongation') ylim([-0.1 1.1]) % set up the arch deformation yyaxis right plot(results_norm.angles_filt.cal_mt1.DP(12,:),'linewidth',3) ylabel(sprintf('Arch dorsiflexion [%s]',0176)) xlabel('% Stance') legend(ligs) legend('Long Plantar','Short Plantar','Calcaneonavicular- Inferior','Calcaneonavicular- Medial','Plantar Fascia','Arch Deformation')