First commit

Committing all code files

Analysis/Analysis_multtrialNN.m

@@ -0,0 +1,162 @@
+clear all
+close all
+clc
+
+%% Add to path
+addpath(genpath('C:\Users\annek\Documents\MATLAB\MuscleNN\musclenn\PaperCode')) %Change this to the folder where you downloaded this code
+
+%% Settings
+
+bird_names = {'Bl3', 'BL4', 'Or3', 'Ye3', 'Pu1'}; %Birds in order of 1 to 5
+muscle_names = {'LG', 'DF'};
+
+%Write here the location of the data
+folder_base = 'C:\Users\annek\Documents\MATLAB\MuscleModel\MuscleData\Guinnea Fowls Share\';
+
+bird_data = readmatrix([folder_base 'MuscleMorphologyData']);
+
+model_params = {'network_large'};% %use 'normal' to compare to original Hill parameters and 'optresult_150523_r1.mat', 'optresult_150523_r12.mat' to compare to optimized Hill models
+model_names = {'NN'};
+
+for k = 1:length(bird_names)
+    trial_names = dir([folder_base filesep bird_names{k}]);
+    % Find some indices for plotting
+    if strcmpi(bird_names{k}, 'or3')
+        clear trial_name
+        for iT = 1:length(trial_names)
+            trial_name{iT} = trial_names(iT).name;
+        end
+        Or3_ind = find(~contains(trial_name, 'r08') == 0);
+    elseif strcmpi(bird_names{k}, 'ye3')
+        clear trial_name
+        trial_names = dir([folder_base filesep bird_names{k}]); 
+        for iT = 1:length(trial_names)
+            trial_name{iT} = trial_names(iT).name;
+        end
+        Ye3_ind = find(~contains(trial_name, 'd2_t11') == 0);
+    elseif strcmpi(bird_names{k}, 'pu1')
+        clear trial_name
+        trial_names = dir([folder_base filesep bird_names{k}]); 
+        for iT = 1:length(trial_names)
+            trial_name{iT} = trial_names(iT).name;
+        end
+        Pu1_ind = find(~contains(trial_name, 'r08_3p8') == 0);
+    end
+
+    for l = 1:length(muscle_names)
+        for i = 3:length(trial_names)
+            if ~strcmp(trial_names(i).name(end-2:end), 'mat')
+                continue
+            end 
+            if strcmpi(bird_names{k}, 'ye3') && strcmpi(muscle_names{l}, 'df') %do not use day 1 for DF of ye3
+                if strcmpi(trial_names(i).name(1:5), 'Ye3d1')
+                    continue
+                end
+            end
+            if strcmpi(bird_names{k}, 'ye3') && strcmpi(muscle_names{l}, 'lg') %do not use LG of ye3
+                continue
+            end
+            if strcmpi(bird_names{k}, 'pu1') && strcmpi(muscle_names{l}, 'df') %do not use DF of pu1
+                continue
+            end
+            %% Get forces, calculate RMSEs and correlations
+            [time(i,k,l).t, data(i,k,l).l_ce, data(i,k,l).v_ce, data(i,k,l).EMG, force(i,1,k,l).meas] = loadDataFile(bird_names{k}, muscle_names{l}, trial_names(i).name,folder_base); 
+            for j = 1:length(model_names)
+                [force(i,j,k,l).est,f_max(i,j,k,l)] = getForce(data(i,k,l).l_ce, data(i,k,l).v_ce, data(i,k,l).EMG, model_names{j}, model_params{j}, muscle_names{l}, bird_names{k}, bird_data);
+
+                RMSE(i,j,k,l) = rmse(force(i,j,k,l).est(:), force(i,1,k,l).meas(:)/f_max(i,j,k,l));
+                corrP(i,j,k,l) = corr(force(i,j,k,l).est(:), force(i,1,k,l).meas(:)/f_max(i,j,k,l));
+                [r2(i,j,k,l), RMSE1(i,j,k,l)] = rsquare(force(i,1,k,l).meas(:)/f_max(i,j,k,l),force(i,j,k,l).est(:));
+            end
+        end
+    end
+end
+
+%% Calculate mean RMSE and correlations, not including training data
+for j = 1:length(model_names)
+    for k = 1:length(bird_names)
+        for l = 1:length(muscle_names)
+%             ind_table = (k-1)*length(muscle_names)+l;
+            RMSE_now = RMSE(:,j,k,l);
+            corrP_now = corrP(:, j,k,l);
+            RMSE1_now = RMSE1(:,j,k,l);
+            R2_now = r2(:, j,k,l);
+
+
+            % identify and remove zeros
+            ind_rem = find(abs(RMSE_now)<1e-6);
+
+            RMSE_now(ind_rem) = [];
+            RMSE_avg(j,k,l) = mean(RMSE_now);
+
+            corrP_now(ind_rem) = [];
+            corrP_avg(j,k,l) = mean(corrP_now);
+
+            RMSE1_now(ind_rem) = [];
+            RMSE1_avg(j,k,l) = mean(RMSE1_now);
+
+            R2_now(ind_rem) = [];
+            R2_avg(j,k,l) = mean(R2_now);
+        end
+    end
+end
+
+%% Plots
+trial_names = {'Or3d1_r08_3p5_7cm_Cal', 'Or3d1_r08_3p5_7cm_Cal', 'Pu1d1_r08_3p8_7cm_Cal', 'ye3d2_t11_3p5_7cm_Cal'};
+
+i1 = [Or3_ind Or3_ind Pu1_ind Ye3_ind];
+k1 = [3 3 5 4];
+l1 = [1 2 1 2]; 
+
+figure
+for i = 1: length(trial_names)
+    subplot(length(trial_names),1,i)
+    plot(time(i1(i),k1(i),l1(i)).t, force(i1(i),1,k1(i),l1(i)).meas, 'k')
+    hold on
+    for j = 1:length(model_names)
+        plot(time(i1(i),k1(i),l1(i)).t, force(i1(i),j,k1(i),l1(i)).est*f_max(i1(i),j,k1(i),l1(i)))
+    end
+    xlim([10 12])
+end
+legend({'Measurements', model_names{:}})
+
+%% Check if force-length and force-velocity relationships are followed
+model_params = {'network_large', 'normal'};
+model_names = {'NN', 'Hill'};
+act = 0.2:0.2:1;
+for i = 1:length(act)
+    l_ce1 = 0.6:0.01:1.4;
+    l_ce1 = l_ce1';
+    for j = 1:length(model_names)
+            flce(:,j,i) = getForce(l_ce1, zeros(size(l_ce1)), act(i)+zeros(size(l_ce1)), model_names{j}, model_params{j}, muscle_names{1});
+    end
+
+    v_ce1 = -10:0.1:10;
+    v_ce1 = v_ce1';
+    for j = 1:length(model_names)
+            gvce(:,j,i) = getForce(ones(size(v_ce1)), v_ce1, act(i)+zeros(size(v_ce1)), model_names{j}, model_params{j}, muscle_names{1});
+    end
+end
+
+figure
+subplot(1,2,1)
+hold on
+for i = 1:length(act)
+    for j = 1:length(model_names)
+        plot(l_ce1, flce(:,j,i))
+    end
+end
+xlabel('Normalized fibre length')
+ylabel('Normalized muscle force')
+title('Force-Length Relationship')
+
+subplot(1,2,2)
+hold on
+for i = 1:length(act)
+    for j = 1:length(model_names)
+        plot(v_ce1, gvce(:,j,i))
+    end
+end
+legend('Hill-type model', 'Neural Network')
+xlabel('Normalized fibre velocity')
+title('Force-Velocity Relationship')

Analysis/Analysis_singletrialNNvsHill.m

@@ -0,0 +1,116 @@
+clear all
+close all
+clc
+
+%% Add to path
+addpath(genpath('C:\Users\annek\Documents\MATLAB\MuscleNN\musclenn\PaperCode')) %Change this to the folder where you downloaded this code
+
+%% Settings
+
+bird_names = {'Bl3', 'BL4'};
+muscle_names = {'LG', 'DF'};
+
+%Write here the location of the data
+folder_base = 'C:\Users\annek\Documents\MATLAB\MuscleModel\MuscleData\Guinnea Fowls\'; %Change this to the folder where the guinnea fowl data is available
+
+bird_data = xlsread([folder_base 'MuscleMorphologyData']);
+
+model_params = {'optresult_150523_r1.mat', 'optresult_150523_r12.mat', 'network_r01' 'network_r12'};
+model_names = {'Hill', 'Hill', 'NN', 'NN'};
+
+for k = 1:length(bird_names)
+    trial_names = dir([folder_base filesep bird_names{k}]);
+    % Find some indices for average calculations and plotting
+    if strcmpi(bird_names{k}, 'bl3')
+        clear trial_name
+        trial_names = dir([folder_base filesep bird_names{k}]); 
+        for iT = 1:length(trial_names)
+            trial_name{iT} = trial_names(iT).name;
+        end
+        r01_ind = find(~contains(trial_name, 'r01') == 0);
+        r10_ind = find(~contains(trial_name, 'r10') == 0);
+        r12_ind = find(~contains(trial_name, 'r12') == 0);
+    elseif strcmpi(bird_names{k}, 'bl4')
+        clear trial_name
+        trial_names = dir([folder_base filesep bird_names{k}]); 
+        for iT = 1:length(trial_names)
+            trial_name{iT} = trial_names(iT).name;
+        end
+        r09_ind = find(~contains(trial_name, 'r09') == 0);
+    end
+
+    for l = 1:length(muscle_names)
+        for i = 3:length(trial_names)
+            if ~strcmp(trial_names(i).name(end-2:end), 'mat')
+                continue
+            end 
+
+            %% Get forces, calculate RMSEs and correlations
+            [time(i,k,l).t, data(i,k,l).l_ce, data(i,k,l).v_ce, data(i,k,l).EMG, force(i,1,k,l).meas] = loadDataFile(bird_names{k}, muscle_names{l}, trial_names(i).name,folder_base); 
+            for j = 1:length(model_names)
+                [force(i,j,k,l).est,f_max(i,j,k,l)] = getForce(data(i,k,l).l_ce, data(i,k,l).v_ce, data(i,k,l).EMG, model_names{j}, model_params{j}, muscle_names{l}, bird_names{k}, bird_data);
+
+                RMSE(i,j,k,l) = rmse(force(i,j,k,l).est(:), force(i,1,k,l).meas(:)/f_max(i,j,k,l));
+                corrP(i,j,k,l) = corr(force(i,j,k,l).est(:), force(i,1,k,l).meas(:)/f_max(i,j,k,l));
+                [r2(i,j,k,l), RMSE1(i,j,k,l)] = rsquare(force(i,1,k,l).meas(:)/f_max(i,j,k,l),force(i,j,k,l).est(:));
+            end
+        end
+    end
+end
+
+%% Calculate mean RMSE and correlations, not including training data
+for j = 1:length(model_names)
+    for k = 1:length(bird_names)
+        for l = 1:length(muscle_names)
+            ind_table = (k-1)*length(muscle_names)+l;
+            RMSE_now = RMSE(:,j,k,l);
+            corrP_now = corrP(:, j,k,l);
+            RMSE1_now = RMSE1(:,j,k,l);
+            R2_now = r2(:, j,k,l);
+
+            % identify and remove zeros
+            ind_rem = find(abs(RMSE_now)<1e-6);
+            %remove training data
+            if strcmpi(bird_names{k}, 'bl3')
+                cur_modname = model_params{j};
+                r1_loc = strfind(cur_modname, 'r1');
+                if strcmpi(cur_modname(r1_loc+2), '2')
+                    ind_rem = [ind_rem; r12_ind];
+                else
+                    ind_rem = [ind_rem; r01_ind];
+                end
+            end
+
+            RMSE_now(ind_rem) = [];
+            RMSE_avg(ind_table,j) = mean(RMSE_now);
+
+            corrP_now(ind_rem) = [];
+            corrP_avg(ind_table,j) = mean(corrP_now);
+
+            RMSE1_now(ind_rem) = [];
+            RMSE1_avg(ind_table,j) = mean(RMSE1_now);
+
+            R2_now(ind_rem) = [];
+            R2_avg(ind_table,j) = mean(R2_now);
+        end
+    end
+end
+
+
+%% Plots
+trial_names = {'Bl3d2_r10_3p8_7cm_Cal', 'Bl3d2_r10_3p8_7cm_Cal', 'BL4d2_r09_3p8_7cm_Cal', 'BL4d2_r09_3p8_7cm_Cal'};
+i1 = [r10_ind r10_ind r09_ind r09_ind];
+k1 = [1 1 2 2];
+l1 = [1 2 1 2];
+
+
+figure
+for i = 1: length(trial_names)
+    subplot(length(trial_names),1,i)
+    plot(time(i1(i),k1(i),l1(i)).t, force(i1(i),1,k1(i),l1(i)).meas, 'k')
+    hold on
+    for j = 1:length(model_names)
+        plot(time(i1(i),k1(i),l1(i)).t, force(i1(i),j,k1(i),l1(i)).est*f_max(i1(i),j,k1(i), l1(i)))
+    end
+    xlim([7.4 8.9])
+end

Analysis/getForce.m

@@ -0,0 +1,110 @@
+function [Force, Force_iso] = getForce(l_ce, v_ce, EMG, model_name, model_params, muscle_name, bird_name, bird_data)
+
+if nargin > 6
+    %% Load parameters from excel
+    ind_row = findBird(bird_name);
+
+    warning('The columns of the xlsx file are hard coded, please make sure that they match your version of MuscleMorphologyData.xlsx')
+
+    if strcmpi(muscle_name, 'lg')
+%         musvar.penn_ang = bird_data(ind_row,7);
+        musvar.PCSA = bird_data(ind_row,8)/1000/1000; %converted to m
+%         musvar.mmass = bird_data(ind_row,5)/1000;
+        musvar.l_opt = bird_data(ind_row,6)/1000;
+    elseif strcmpi(muscle_name, 'df')
+%         musvar.penn_ang = bird_data(ind_row,13);
+        musvar.PCSA = bird_data(ind_row,14)/1000/1000;
+%         musvar.mmass = bird_data(ind_row,11)/1000;
+        musvar.l_opt = bird_data(ind_row,12)/1000;
+    else
+        error('Incorrect muscle name')
+    end
+
+    sigma = 3.6e5; % Max muscle stress, from: S M Cox, K L Easton, M Cromie Lear, R L Marsh, S L Delp, J Rubenson, The Interaction of Compliance and Activation on the Force-Length Operating Range and Force Generating Capacity of Skeletal Muscle: A Computational Study using a Guinea Fowl Musculoskeletal Model, Integrative Organismal Biology, Volume 1, Issue 1, 2019, obz022, https://doi.org/10.1093/iob/obz022
+    musvar.f_max = sigma*musvar.PCSA;
+else
+    musvar.f_max = 1;
+    musvar.l_opt = 1;
+end
+
+if strcmpi(model_name, 'hill')
+    if strcmpi(model_params(end-3:end), '.mat')
+
+        load(model_params) %Change file name to change model parameters
+        %assign optimized variables
+        modelvar.v_max = optvar.Position(1);
+        modelvar.Arel = optvar.Position(2);
+        modelvar.gmax = optvar.Position(3);
+        modelvar.W = optvar.Position(4);
+        modelvar.kPEE = optvar.Position(5);
+        modelvar.PEEslack = optvar.Position(6);
+    else
+        modelvar.v_max = 10;
+        modelvar.PEEslack = 1.2;
+        modelvar.gmax = 1.5;
+        modelvar.kPEE = 1/musvar.l_opt^2;
+        modelvar.Arel = 0.25;
+        modelvar.W = 0.4;
+    end
+
+    %dependent parameters         
+    modelvar.c3 = modelvar.v_max*modelvar.Arel*(modelvar.gmax - 1.)/(modelvar.Arel + 1);
+
+    %% Calculate Hill forces
+    if abs(std(l_ce)) < 1e-4 %means to check force-velocity relationship
+        Force = findMuscleForce(modelvar, musvar, ones(size(v_ce)), v_ce, EMG, 'gvce');
+    elseif abs(std(v_ce)) < 1e-4 %means to check force-length relationship including PEE
+        Force = findMuscleForce(modelvar, musvar, l_ce, zeros(size(l_ce)), EMG, 'flcePEE');
+    else
+        Force = findMuscleForce(modelvar, musvar, l_ce, v_ce, EMG);
+    end
+elseif strcmpi(model_name, 'nn')
+    load(model_params) %Change file name to change model parameters
+
+    Mdl = NN.Mdl;
+    data_mean = NN.data_mean;
+    data_std = NN.data_std;
+
+    if abs(std(l_ce)) < 1e-4 %means to check force-velocity relationship
+        lce_norm = doDataNormalization(ones(size(v_ce)), data_mean.lce,data_std.lce);
+        vce_norm = doDataNormalization(v_ce, data_mean.vce,data_std.vce);
+        EMG_norm = doDataNormalization(EMG, data_mean.EMG,data_std.EMG);
+    elseif abs(std(v_ce)) < 1e-4 %means to check force-length relationship
+        lce_norm = doDataNormalization(l_ce, data_mean.lce,data_std.lce);
+        vce_norm = doDataNormalization(zeros(size(l_ce)), data_mean.vce,data_std.vce);
+        EMG_norm = doDataNormalization(EMG, data_mean.EMG,data_std.EMG);
+    elseif abs(std(EMG)) < 1e-4 %means to check PEE relationship
+        lce_norm = doDataNormalization(l_ce, data_mean.lce,data_std.lce);
+        vce_norm = doDataNormalization(zeros(size(l_ce)), data_mean.vce,data_std.vce);
+        EMG_norm = doDataNormalization(EMG, data_mean.EMG,data_std.EMG);
+    else      
+        lce_norm = doDataNormalization(l_ce, data_mean.lce,data_std.lce);
+        vce_norm = doDataNormalization(v_ce, data_mean.vce,data_std.vce);
+        EMG_norm = doDataNormalization(EMG, data_mean.EMG,data_std.EMG);
+    end
+
+    X = [lce_norm vce_norm EMG_norm];
+    %% Calculate NN forces
+    Force_norm = predict(Mdl, X);
+    Force = doDataDenormalization(Force_norm, data_mean.force, data_std.force);
+elseif strcmpi(model_name, 'nn_emg')
+    load(model_params) %Change file name to change model parameters
+
+    Mdl = NN.Mdl;
+    data_mean = NN.data_mean;
+    data_std = NN.data_std;
+
+    EMG_norm = doDataNormalization(EMG, data_mean.EMG,data_std.EMG);
+
+    X = EMG_norm;
+    %% Calculate NN forces
+    Force_norm = predict(Mdl, X);
+    Force = doDataDenormalization(Force_norm, data_mean.force, data_std.force);
+else
+    error('incorrect model_name, should be Hill or NN')
+end
+
+if nargout > 1
+    Force_iso = musvar.f_max;
+end
+