Evaluating parameter value estimation accuracy - simulated data

The analysis of this page corresponds to the section with the same title (i.e., ‘Evaluating parameter value estimation accuracy - simulated data’).

There are three steps:

Step 1: 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 2: 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 vPar in ['TM', 'IM']:
        # Load parameters
        parFile = os.path.join(dataDir, mus, 'parameters', f'{mus}_{vPar}.pkl')
        muspar = pickle.load(open(parFile, 'rb'))[0]
        for exp in ['QR','SR','ISOM']:
            dataDirExp = os.path.join(dataDir, mus, 'dataExp', exp)
            dataDirSim = os.path.join(dataDir, mus, 'simsExp', vPar, exp, '')

            # 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 = simulate_file(filepath, muspar)
                fname = os.path.basename(filepath)[:-6]
                outName = f"{fname}{vPar}.csv"
                # df.to_csv(os.path.join(dataDirSim, outName), index=False)
                print(f"Saved: {outName}")

Step 3: Compute variables for manuscript

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

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

# Custom functions
import hillmodel, stats

# %% Extract parameters
muscles = ['GMs1','GMs2','GMs3']
param_keys = ['a', 'b', 'kpee', 'ksee', 'fmax', 'lce_opt', 'lpee0', 'lsee0', 'tact', 'tdeact']

orParms, tmParms, imParms = [], [], []
for mus in muscles:   
    orPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl'), 'rb'))
    tmPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))[0]
    imPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_IM.pkl'), 'rb'))[0]

    parms = [orPar[k] for k in param_keys]
    orParms.append(parms)
    
    parms = [tmPar[k] for k in param_keys]
    tmParms.append(parms)
    
    parms = [imPar[k] for k in param_keys]
    imParms.append(parms)
    
# Store in pandas dataframe
df_or = pd.DataFrame(list(zip(*orParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with actual/real parameter values
df_tm = pd.DataFrame(list(zip(*tmParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with estimated parameter values using the traditional method
df_im = pd.DataFrame(list(zip(*imParms)), columns=['GM1', 'GM2', 'GM3'], index=param_keys) # df with estimated parameter values using the improved method

# %% Compute percentage differences from real/actual and estimated parameter values
# Positive: Estimated is that % higher than Real
# Negative: Estimated is that % lower than Real
df_tm_p = stats.pdiff(df_tm, df_or)
df_im_p = stats.pdiff(df_im, df_or)

# %% Variables related to overestimation of length changes due to QR
deltaLsee_tm   = np.full((len(muscles),10), np.nan)
deltaLsee_real = np.full((len(muscles),10), np.nan)
for iMus,mus in enumerate(muscles):   
    orPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl'), 'rb'))
    tmPar,dataQRtm = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))[0:2]

    fseeQRpre = dataQRtm['fseeQRpre']
    fseeQRpst = dataQRtm['fseeQRpst']
    lseeQRpre = (fseeQRpre/orPar['ksee'])**0.5 + orPar['lsee0']
    lseeQRpst = (fseeQRpst/orPar['ksee'])**0.5 + orPar['lsee0']
    
    # Compute change in SEE length due to QR
    deltaLsee_tm[iMus] = dataQRtm['lseeQRpre']-dataQRtm['lseeQRpst']  # [m] estimated SEE length change due to QR
    deltaLsee_real[iMus] = lseeQRpre-lseeQRpst  # [m] real/actual SEE length change due to QR

# Compute percentage difference from real to estimated SEE length change due to QR
# Positive: TM is that % higher than Real
# Negative: TM is that % lower than Real
p_deltaLsee = stats.pdiff(deltaLsee_tm, deltaLsee_real)  # [%] percentage difference 
p_deltaLseeAvg = np.mean(p_deltaLsee,1)  # [%] average per muscle
p_deltaLseeStd = np.std(p_deltaLsee,1)  # [%] std per muscle

# Compute difference between real and estimated CE length change due to QR
deltaLce = deltaLsee_real-deltaLsee_tm  # [m]
deltaLceAvg = np.mean(deltaLce,1)  # [m] average per muscle
deltaLceStd = np.std(deltaLce,1)  # [m] std per muscle

# %% 3.1.1: Traditional vs. Improved method
# Extract all contraction dynamics parameters except ksee
p_TM_con = df_tm_p.loc[['a','b', 'kpee', 'fmax', 'lce_opt', 'lpee0', 'lsee0']]

# Avg. difference: from real to estimated SEE length change
p_TM_deltaLsee_avgall = f'{p_deltaLseeAvg.mean():0.0f}' #  [%] avg. percentage difference
# Avg. difference over all muscles: from real to estimated SEE stiffness
p_TM_ksee_avg = f'{df_tm_p.loc['ksee'].mean():0.0f}' #  [%] avg. percentage difference
# Avg. difference over all muscles: from real to estimated tact
p_TM_tact_avg = f'{df_tm_p.loc['tact'].mean():0.0f}'  # [%] avg. percentage difference
# Max. difference over all muscles: from real to estimated tact
p_TM_con_max = f'{p_TM_con.abs().max().max():0.1f}'  # [%] max. percentage difference
# Max. difference over all muscles: from real to estimated tact
p_IM_max = f'{df_im_p.abs().max().max():0.0f}'  # [%] max. percentage difference

# %% 3.1.2: Estimation of SEE stiffness. 
# Difference between SEE elongation of estimated parameter value and real/actual value
# Positive: estimated value is higher than real/actual value
d_esee = (df_tm.loc['fmax']/df_tm.loc['ksee'])**0.5 - (df_or.loc['fmax']/df_or.loc['ksee'])**0.5 # [m]

# Correlation coefficients
r_see = np.full(len(muscles), np.nan)
r_pee = np.full(len(muscles), np.nan)
r_mtc = np.full(len(muscles), np.nan)
for iMus,mus in enumerate(muscles):   
    parFile = os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl')
    muspar, dataQR, dataSR, dataACTout = pickle.load(open(parFile, 'rb'))
    
    # SEE
    lseeData = np.hstack((dataQR['lseeQRpre'],dataQR['lseeQRpst']))
    fseeData = np.hstack((dataQR['fseeQRpre'],dataQR['fseeQRpst']))
    fseeMdl = hillmodel.LEE2Force(lseeData,0,muspar)[0]
    r_see[iMus] = np.corrcoef(fseeData, fseeMdl)[0,1]
    
    # PEE
    fpeeMdl = hillmodel.LEE2Force(0,dataQR['lpeeQR'],muspar)[1]
    r_pee[iMus] = np.corrcoef(dataQR['fpeeQR'], fpeeMdl)[0,1]
    
    # MTC
    fseeMdl = hillmodel.ForceEQ(dataQR['lmtcQRpre'],1,muspar)[0]
    r_mtc[iMus] = np.corrcoef(dataQR['fseeQRpre'], fseeMdl)[0,1]
r_all = np.hstack((r_see,r_pee,r_mtc))

# Percentage difference from real to estimated parameter
p_TM_ksee = [f'{v:0.0f}' for v in df_tm_p.loc['ksee']]  # [%] per muscle
# Avg&Std per muscle: differences from real to estimated parameters 
p_TM_deltaLsee_avgmus = [f'{v:0.0f}' for v in p_deltaLseeAvg]  # [%]
p_TM_deltaLsee_stdmus = [f'{v:0.0f}' for v in p_deltaLseeStd]  # [%]
# Avg&Std per muscle: CE shortening during QR
p_TM_deltaLce_avgmus = [f'{v*1e6:0.0f}' for v in deltaLceAvg]  # [um] avg. per muscle
p_TM_deltaLce_stdmus = [f'{v*1e6:0.0f}' for v in deltaLceStd]  # [um] std per muscle
# Difference from real to estimate in SEE elongation @ fmax
d_TM_esee = [f'{v*1e3:0.2f}' for  v in d_esee] #  [mm] 
# Minimum correlation coefficient
r_TM_min = f'{np.min(np.floor(r_all*100)/100):0.2f}'  # []

# %% 3.1.2: Estimation of PEE parameter values
nPEEdata = np.full(len(muscles), np.nan)
for iMus,mus in enumerate(muscles):   
    muspar, dataQR, dataSR, dataACT = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM.pkl'), 'rb'))
    lpeeQR = dataQR['lpeeQR']  # [m] PEE length of data
    fpeeQR = dataQR['fpeeQR']  # [N] PEE force of data
    fpeeModel = hillmodel.LEE2Force(0,lpeeQR,muspar)[1]  # [N] PEE force estimated
    
    # Compute how many data points available
    nPEEdata[iMus] = np.sum(fpeeModel>0)

# Dislay percentage differences from real to estimated parameters
p_TM_kpee = [f'{v:0.1f}' for v in df_tm_p.loc['kpee']]  # [%] percentage difference
p_TM_lpee0 = [f'{v:0.1f}' for v in df_tm_p.loc['lpee0']]  # [%] percentage difference

# Display amount of datapoints available
n_PEE_data = [f'{v:0.0f}' for v in nPEEdata]  # [ ]

# %% 3.1.2: Estimation of Fcemax, Lceopt, Lsee0
# Dislay percentage differences from real to estimated parameters
p_TM_fl_max = f'{df_tm_p.loc[['fmax','lce_opt','lsee0'],:].abs().max().max():0.1f}'  # [%] maximal percentage difference

# %% 3.1.2: Estimation of CE force-velocity parameter values
# Dislay percentage differences from real to estimated parameters
p_TM_a_max = f'{df_tm_p.loc['a',:].abs().max():0.1f}'  # [%] maximal percentage difference
p_TM_b_max = f'{df_tm_p.loc['b',:].abs().max():0.1f}'  # [%] maximal percentage difference

# %% 3.1.2: Estimation of excitation dynamics parameter values
p_tact_QR = np.full(len(muscles), np.nan)
p_tdeact_QR = np.full(len(muscles), np.nan)
p_tact_SR = np.full(len(muscles), np.nan)
p_tdeact_SR = np.full(len(muscles), np.nan)
for iMus, mus in enumerate(muscles):   
    orPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl'), 'rb'))
    qrPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM_QR.pkl'), 'rb'))[0]
    srPar = pickle.load(open(os.path.join(dataDir,mus,'parameters',mus+'_TM_SR.pkl'), 'rb'))[0]
    
    # Compute percentage difference from real to estimated parameter
    # Positive: estimated is that % higher than Real
    # Negative: estimated is that % lower than Real
    p_tact_QR[iMus] = stats.pdiff(qrPar['tact'],orPar['tact'])
    p_tact_SR[iMus] = stats.pdiff(srPar['tact'],orPar['tact'])
    p_tdeact_QR[iMus] = stats.pdiff(qrPar['tdeact'],orPar['tdeact'])
    p_tdeact_SR[iMus] = stats.pdiff(srPar['tdeact'],orPar['tdeact'])

# Dislay percentage differences from real to estimated parameters by using ISOM data
p_TM_tact = [f'{v:0.0f}' for v in df_tm_p.loc['tact']]  # [%] percentage difference
p_TM_tdeact = [f'{v:0.1f}' for v in df_tm_p.loc['tdeact']]  # [%] percentage difference

# Display percentage difference from real to estimated parameters by using
# either QR or SR data
p_tact_QR_avg = f'{p_tact_QR.mean():0.0f}'  #  [%] percentage difference when using QR data
p_tact_SR_avg = f'{p_tact_SR.mean():0.0f}'  #  [%] percentage difference when using SR data
p_tdeact_QR_avg = f'{p_tdeact_QR.mean():0.0f}'  #  [%] percentage difference when using QR data
p_tdeact_SR_avg = f'{p_tdeact_SR.mean():0.0f}'  #  [%] percentage difference when using SR data
# For tdeact: differences are similar so show avg of the two
p_tdeact_QRSR = np.concatenate([p_tdeact_QR, p_tdeact_SR])
p_tdeact_QRSR_avg = f'{p_tdeact_QRSR.mean():0.0f}'  #  [%] percentage difference when using QR/SR data