Sensitivity analysis

The analysis of this page corresponds to the section with the same title (i.e., ‘Sensitivity analysis’). There are two sensivitiy analysis performed:

Monte Carlo simulations
Interdependency of parameter values

Monte Carlo simulations

Three steps are taken for the Monte Carlo simulation:

Step 1A: Estimate muscle property values.
Step 1B: Simulate protocols to predict SEE force over time.
Step 1C: Compute variables for manuscript.

Step 1A: Estimate muscle property values


# %% Load packages
import os, pickle 
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent
dataDir = baseDir / 'data'

# Load custom functions
from functions.mp_estimator import loaddata, estimate

# %% Muscle & folders
for mus in ['GMs1', 'GMs2', 'GMs3']: 
    dataDirMus = os.path.join(dataDir,mus,'dataExp','')
    parFile = os.path.join(dataDir,mus,'parameters',mus)
    
    # %% Settings and load data
    optsQR = {
        'dataDir':          dataDirMus+'QR',
        'iCols':            (0,1,3),
        'idxQRmin':         'auto',
        'idxQRpre':         'auto',
        'idxQRpst':         'auto',
        'nQRsamp':          'auto',
        'dispFig':          False,
        }
    optsSR = {
        'dataDir':          dataDirMus+'SR',
        'iCols':            (0,1,3),
        'idxSRcon':         'auto',
        'nSRsamp':          'auto',
        'dispFig':          False,
        }
    optsISOM = {
        'dataDir':          dataDirMus+'SR',
        'iCols':            (0,1,3,2),
        'durStimOffset':    0.1,
        'dispFig':          False,
        }
    dataQR,idxQRmin,idxQRpre,idxQRpst = loaddata.qr(optsQR)
    dataSR,idxSRcon = loaddata.sr(optsSR)
    dataISOM, idxSEL = loaddata.isom(optsISOM)
    
    # %% Traditional method
    muspar = loaddata.get_muspar()
    estpar,dataQRout = estimate.fl(dataQR,muspar)
    estpar,dataSRout = estimate.fv(dataSR,muspar,estpar)
    estpar,dataACTout = estimate.act(dataISOM,muspar,estpar,do_print=True)   
    pickle.dump([{**muspar,**estpar}, dataQRout, dataSRout, dataACTout], open(parFile+'_TM_SR.pkl', 'wb'))
    
    # %% Improved method
    muspar = loaddata.get_muspar()
    estparr,dataQRout,dataSRout = estimate.im(dataQR,dataSR,muspar,do_print=True)
    estparr,dataACTout = estimate.act(dataISOM,muspar,estparr,do_print=True)
    pickle.dump([{**muspar,**estparr}, dataQRout, dataSRout, dataACTout], open(parFile+'_IM.pkl', 'wb'))

Step 1B: Simulate protocols to predict SEE force over time

# %% Load packages & set directories
import pickle, os, glob
from pathlib import Path
from functions.helpers import simulate_file

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent
dataDir = baseDir / 'data'

# %% Simulate to predict SEE force based on experimental data
muscles = ['GMs1', 'GMs2', 'GMs3']

for mus in muscles:
    for exp in ['QR','SR','ISOM']:
        for iMC in range(1,51):
            parFile = os.path.join(dataDir,mus,'parameters','mc',f'{mus}_MC{iMC:02d}'+'.pkl')
            muspar = pickle.load(open(parFile, 'rb'))[0]
            
            dataDirExp = os.path.join(dataDir,mus,'dataMC',exp,f'{iMC:02d}','')
            dataDirSim = os.path.join(dataDir,mus,'simsMC',exp,f'{iMC:02d}','')

            # Clean old files
            for f in glob.glob(os.path.join(dataDirSim, '*')):
                os.remove(f)
            
            # Process each file
            for filepath in glob.glob(os.path.join(dataDirExp, '*')):
                df_sim, df_exp = simulate_file(filepath, muspar)
                fname = os.path.basename(filepath)[:-4]
                outName = f"{fname}.csv"
                df_sim.to_csv(os.path.join(dataDirSim, outName), index=False)
                print(f"Saved: {outName}")

Step 1C: Compute variables for manuscript

# %% Load packages
import os, pickle, sys
import pandas as pd
from pathlib import Path

# Set directories
cwd = Path.cwd()
baseDir = cwd.parent
dataDir = baseDir / 'data'
funcDir = baseDir / 'analysis' / 'functions'
sys.path.append(str(funcDir))

# Custom functions
import stats

# %% 3.2.1: Monte Carlo simulations
muscles = ['GMs1', 'GMs2', 'GMs3']
# Parameters to track
params = ['a', 'b', 'fmax', 'kpee', 'ksee', 'lce_opt', 'lpee0', 'lsee0', 'tact', 'tdeact']

# Store Monte-Carlo results
results = []
for mus in muscles:
    par_file = os.path.join(dataDir, mus, 'parameters', f'{mus}_OR.pkl')
    or_par = pickle.load(open(par_file, 'rb'))
    eseeMaxOR = (or_par['fmax']/or_par['ksee'])**0.5
        
    
    for iMC in range(1, 51):
        mc_file = os.path.join(dataDir, mus, 'parameters','mc', f'{mus}_MC{iMC:02d}.pkl')
        mc_par = pickle.load(open(mc_file, 'rb'))[0]
        eseeMaxMC = (mc_par['fmax']/mc_par['ksee'])**0.5
        
        entry = {'muscle': mus, 'MC': iMC, 'd_esee': (eseeMaxMC-eseeMaxOR)}
        for param in params:
            entry[param] = stats.pdiff(mc_par[param], or_par[param])
        
        results.append(entry)

# Convert to pandas DataFrame
df = pd.DataFrame(results)

# Compute summary statistics
df_mean = df.groupby('muscle').mean().T
df_std = df.groupby('muscle').std().T

# Extract all contraction dynamics parameters except ksee
p_MC_con_avg = df_mean.loc[['a','b', 'kpee', 'fmax', 'lce_opt', 'lpee0', 'lsee0']]
p_MC_con_std = df_std.loc[['a', 'b', 'fmax', 'kpee', 'lce_opt', 'lpee0', 'lsee0']]

# Percentage difference over all muscles: from real to estimated ksee
p_MC_ksee_avg = f'{df_mean.loc['ksee'].mean():0.0f}'  # [%] percentage difference
# Absolute difference over all muscles of SEE elongation @ Fcemax
d_MC_esee_avg = f'{df_mean.loc['d_esee'].mean()*1e3:0.1f}'  # [mm]
# Max. avg. percentage difference over all muscles: from real to estimated value of contraction dyn. parms
p_MC_con_maxavg = f'{p_MC_con_avg.abs().max().max():0.1f}' 
# Max. std. percentage difference over all muscles: from real to estimated value of contraction dyn. parms
p_MC_con_maxstd = f'{p_MC_con_std.abs().max().max():0.0f}' 
# Max. std percentage difference over all muscles: from real to estimated value of kpee
p_MC_kpee_maxstd = f'{df_std.loc['kpee'].max():0.0f}'
# Max. avg percentage difference over all muscles: from real to estimated value of tact
p_MC_tact_maxavg = f'{df_mean.loc['tact'].max():0.1f}'
# Avg. std percentage difference over all muscles: from real to estimated value of tact
p_MC_tact_avgstd = f'{df_std.loc['tact'].mean():0.0f}'

Interdependency of parameter values

# %% Load packages
import os, pickle

# Set directories
cwd = os.path.dirname(os.path.abspath(__file__))
baseDir = os.path.join(cwd,'..')
dataDir = os.path.join(baseDir,'data')

# Load custom functions
from functions.mp_estimator import loaddata, estimate_id

#%% Muscle & folders
for iMus, mus in enumerate(['GMs1','GMs2','GMs3']): 
    dataDirMus = os.path.join(dataDir,mus,'dataExp','')
    parDir = os.path.join(dataDir,mus,'parameters','')
    
    #%% Settings and load data
    optsQR = {
        'dataDir':          dataDirMus+'QR',
        'iCols':            (0,1,3),
        'idxQRmin':         'auto',
        'idxQRpre':         'auto',
        'idxQRpst':         'auto',
        'nQRsamp':          'auto',
        'dispFig':          False,
        }
    optsSR = {
        'dataDir':          dataDirMus+'SR',
        'iCols':            (0,1,3),
        'idxSRcon':         'auto',
        'nSRsamp':          'auto',
        'dispFig':          False,
        }
    optsISOM = {
        'dataDir':          dataDirMus+'ISOM',
        'iCols':            (0,1,3,2),
        'durStimOffset':    0.1,
        'dispFig':          False,
        }
    dataQR,idxQRmin,idxQRpre,idxQRpst = loaddata.qr(optsQR)
    dataSR,idxSRcon = loaddata.sr(optsSR)
    dataISOM, idxSEL = loaddata.isom(optsISOM)
    
    #%% Perform sensitivity analyis
    #for sPar in ['a', 'b', 'fmax', 'ksee', 'lce_opt', 'lsee0']:
    for sPar in ['ksee', 'lce_opt', 'lsee0']:
        for idx,fChange in enumerate([0.95, 1.05]):
            muspar = pickle.load(open(parDir+mus+'_IM.pkl', 'rb'))[0]
            print(muspar[sPar])
            muspar[sPar] = muspar[sPar]*fChange
            print(muspar[sPar])
            
            estpar,dataQRout,dataSRout = estimate_id.im(dataQR,dataSR,muspar,{sPar: muspar[sPar]},do_print=True)  
            estpar,dataACTout = estimate_id.act(dataISOM,muspar,estpar,do_print=True)
            filepath = os.path.join(parDir,'interdep',mus+'_'+sPar+f'_{fChange*100:03.0f}.pkl')
            pickle.dump([{**muspar,**estpar}, dataQRout, dataSRout, dataACTout], open(filepath, 'wb'))
            print(estpar[sPar])