Imposed SSCs with constant MTC shortening/lengthening velocity

The analysis of this page corresponds to the sections ‘Imposed stretch-shortening cycles — constant MTC shortening/lengthening velocity’.

Because the analysis is relatively extensive, it is split into separate scripts included on this page:

  1. Predict the maximally attainable AMPO for various combinations of cycle frequency, FTS and MTC length excursion.
  2. Predict the maximally attainable AMPO for one imposed SSC parameter, while optimising the other two.
  3. Analyse the simulation results.

Custom functions used:

Compute AMPO for various combinations of SSC parameters

Code 
"""
This script predicts AMPO for a full grid of constant-velocity SSCs.

For each combination of cycle frequency, FTS and MTC length excursion,
stimulation onset and offset are optimised while MTC shortening and lengthening
velocity remain constant within each phase. The resulting simulations are saved
to `simsCV2` and provide the coarse parameter grid used to initialise the
one-parameter optimisation scripts.
"""

#%% Load packages & set directories
import os, sys, pickle
import numpy as np
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))

import derive_predictions

#%% Load muscle parameters
mus = 'GMe3'
parFile = os.path.join(dataDir,mus,mus+'_IM.pkl')
muspar = pickle.load(open(parFile, 'rb'))[0]

# Now we change lsee0 to 3mm
muspar_lsee0_new = 3e-3
eseerelmax = (muspar['fmax']/muspar['ksee'])**0.5/muspar['lsee0']
muspar['ksee'] = muspar['ksee']*(muspar['lsee0']/muspar_lsee0_new)**2
muspar['lsee0'] = muspar_lsee0_new
lmtc_opt = muspar['lce_opt']+muspar['lsee0']+eseerelmax*muspar['lsee0']
lmtc_avg = lmtc_opt

#%% Compute AMPO across the SSC parameter grid
# problematic... GMe1: cf=5.0 Hz, fts=0.45 & mle=10mm
cf_set = np.arange(0.5,6.1,0.5) # [Hz] n = 12
fts_set = np.arange(0.05,0.96,0.05) # [] n = 19
mle_set = np.arange(2,11.1,1)*1e-3 # [m] n = 10
# thus 2280 per rat!

# spyder 1
# mle_set = mle_set # GMe3
# print(mle_set)

# spyder 2
# mle_set = mle_set[8:] # GMe2 - 3 t/m 7 mm
# fts_set = fts_set[9:]
# # cf_set = cf_set[9:10]
# print(mle_set)
# print(fts_set)
# print(cf_set)

# spyder 3
mle_set = mle_set[-1::-1] # GMe3
print(mle_set)

# spyder 4
# mle_set = mle_set[16:] # GMe2

i = 0
for mle in mle_set:
    for fts in fts_set:
        for cf in cf_set:
            # Perform optimisation
            initialGuess = {}
            initialGuess['stimGuess'] = [-1e-3, 0.5*fts/cf]
            #initialGuess['stimGuess'] = [-1e-3, 0.5*fts/cf+1e-3]
            
            AMPO, y = derive_predictions.opt_stim(None,2,cf,fts,mle,lmtc_avg,muspar,initialGuess)
            time, lmtc, stim, gamma, lcerel, q, lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel = y[0:14]
            
            # Save results            
            filename = mus+f'_cf{cf:0.1f}Hz_fts{fts:0.2f}_mle{mle*1e3:01.1f}mm'
            filepath = os.path.join(dataDir,mus,'simsCV2',filename+'.csv')
            try:
                df = pd.read_csv(filepath)
                data = df.to_numpy().T
                Pmech = np.trapezoid(data[3],data[1])/data[0][-1]
            except:
                Pmech = -np.inf
            
            if AMPO > Pmech:
                data = np.vstack((time,lmtc,stim,fsee,gamma,lcerel,q,fisomrel,fce)).T       
                pd.DataFrame(data).to_csv(filepath,index=False,header=['Time [s]','Lmtc [m]',
                    'STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]','q [ ]',
                    'fisomrel [ ]', 'fce [N]']) 
            i += 1
            print(i)

Compute AMPO for imposed cycle frequency

Code 
"""
This script predicts the maximally attainable AMPO for SSCs with imposed cycle
frequency using constant MTC shortening/lengthening velocity.

For each imposed cycle frequency, the script first loads the predictions for
various combinations of FTS and MTC length excursion. These predictions are
used to make the initial guess and to define local bounds for the optimisation.
Then, FTS, MTC length excursion and stimulation timing are optimised, and the
result is saved to `simsCV`.
"""

#%% Load packages & set directories
import os, sys, pickle
import numpy as np
import pandas as pd
from scipy import integrate
from pathlib import Path

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

import derive_predictions, helpers, stimulation

import matplotlib.pyplot as plt
plt.close('all')

#%% Load muscle parameters
mus = 'GMe1'
parFile = os.path.join(dataDir,mus,mus+'_IM.pkl')
muspar = pickle.load(open(parFile, 'rb'))[0]

# Now we change lsee0 to 3mm
muspar_lsee0_new = 3e-3
eseerelmax = (muspar['fmax']/muspar['ksee'])**0.5/muspar['lsee0']
muspar['ksee'] = muspar['ksee']*(muspar['lsee0']/muspar_lsee0_new)**2
muspar['lsee0'] = muspar_lsee0_new
lmtc_avg = muspar['lce_opt']+(1+eseerelmax)*muspar['lsee0']

#%% Optimise FTS and MTC length excursion for imposed cycle frequency
dataDirSim = os.path.join(dataDir,mus,'simsCV','')

cf_set = np.arange(1, 6.1, 0.5) # [Hz]
for cf in cf_set[0:1]:
    # First read out solution for imposed CF, FTS & MLE
    fts_set = np.arange(0.05, 0.96, 0.05) # []
    mle_set = np.arange(2, 11.1, 1)*1e-3 # [m]
    AMPO_set = helpers.load_sims(cf,fts_set,mle_set,mus,dataDirSim)
        
    # Find maximum
    max_idx = np.nanargmax(AMPO_set)
    iRow, iCol = np.unravel_index(max_idx, AMPO_set.shape)
    
    # Create bounds
    if iRow == 0:
        fts_bounds = (0.1, fts_set[iRow+1])
    elif iRow == fts_set.shape[0]-1:
        fts_bounds = (fts_set[iRow-1], 14)
    else:
        fts_bounds = (fts_set[iRow-1], fts_set[iRow+1])
    if iCol == 0:
        mle_bounds = (1e-3, mle_set[iCol+1])
    elif iCol == mle_set.shape[0]-1:
        mle_bounds = (mle_set[iCol-1], 15e-3)
    else:
        mle_bounds = (mle_set[iCol-1], mle_set[iCol+1])
    
    # Now create initial guess.
    fileName = mus+f'_cf{cf:0.1f}Hz_fts{fts_set[iRow]:0.2f}_mle{mle_set[iCol]*1e3:0.1f}mm'
    df = df = pd.read_csv(dataDirSim+fileName+'.csv')
    data = df.to_numpy()
    time,lmtc,stim,fsee,_,_ = data.T[0:6]
    t_stimOn, t_stimOff = stimulation.get_stim_timing(time,stim)
    t_stimOn = t_stimOn[-1]-time[-1]
    t_stimOff = t_stimOff[-1]
        
    initialGuess = {}
    initialGuess['ftsGuess'] = fts_set[iRow]
    initialGuess['mleGuess'] = mle_set[iCol]
    initialGuess['stimGuess'] = [t_stimOn, t_stimOff]
    
    # Perform optimisation
    AMPO, y, optPar = derive_predictions.opt_ssc_par(2,cf,fts_bounds,mle_bounds,lmtc_avg,muspar,initialGuess)
    time, lmtc, stim, gamma, lcerel, q, lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel = y[0:14]
    
    # Save results
    filename = mus+f'_cf{cf:0.1f}Hz_ftsOpt_mleOpt'
    filepath = os.path.join(dataDir,mus,'simsCV',filename+'.csv')
    try:
        df = pd.read_csv(filepath)
        data = df.to_numpy().T
        Pmech = integrate.trapezoid(data[3],data[1])/data[0][-1]
    except:
        Pmech = 0
    
    if AMPO > Pmech:
        # Save data
        data = np.vstack((time,lmtc,stim,fsee,gamma,lcerel,q,fisomrel,fce)).T
        pd.DataFrame(data).to_csv(filepath,index=False,header=['Time [s]','Lmtc [m]',
            'STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]','q [ ]',
            'fisomrel [ ]', 'fce [N]'])

Compute AMPO for imposed FTS

Code 
"""
This script predicts the maximally attainable AMPO for SSCs with imposed FTS
using constant MTC shortening/lengthening velocity.

For each imposed FTS value, the script first loads the predictions for various
combinations of cycle frequency and MTC length excursion. These predictions are
used to make the initial guess and to define local bounds for the optimisation.
Then, cycle frequency, MTC length excursion and stimulation timing are
optimised, and the result is saved to `simsCV`.
"""

#%% Load packages & set directories
import os, sys, pickle
import numpy as np
import pandas as pd
from scipy import integrate
from pathlib import Path

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

import derive_predictions, helpers, stimulation

import matplotlib.pyplot as plt
plt.close('all')

#%% Load muscle parameters
mus = 'GMe1'
parFile = os.path.join(dataDir,mus,mus+'_IM.pkl')
muspar = pickle.load(open(parFile, 'rb'))[0]

# Now we change lsee0 to 3mm
muspar_lsee0_new = 3e-3
eseerelmax = (muspar['fmax']/muspar['ksee'])**0.5/muspar['lsee0']
muspar['ksee'] = muspar['ksee']*(muspar['lsee0']/muspar_lsee0_new)**2
muspar['lsee0'] = muspar_lsee0_new
lmtc_avg = muspar['lce_opt']+(1+eseerelmax)*muspar['lsee0']

#%% Optimise cycle frequency and MTC length excursion for imposed FTS
dataDirSim = os.path.join(dataDir,mus,'simsCV','')

fts_set = np.arange(0.05, 0.96, 0.05) # []
for fts in fts_set[-1:]:
    # First read out solution for imposed CF, FTS & MLE
    cf_set = np.arange(0.5, 6.1, 0.5) # [Hz]
    mle_set = np.arange(2, 11.1, 1)*1e-3 # [m]
    AMPO_set = helpers.load_sims(cf_set,fts,mle_set,mus,dataDirSim)
        
    # Find maximum
    max_idx = np.nanargmax(AMPO_set)
    iRow, iCol = np.unravel_index(max_idx, AMPO_set.shape)
    
    # Create bounds
    if iRow == 0:
        cf_bounds = (0.1, cf_set[iRow+1])
    elif iRow == cf_set.shape[0]-1:
        cf_bounds = (cf_set[iRow-1], 14)
    else:
        cf_bounds = (cf_set[iRow-1], cf_set[iRow+1])
    if iCol == 0:
        mle_bounds = (1e-3, mle_set[iCol+1])
    elif iCol == mle_set.shape[0]-1:
        mle_bounds = (mle_set[iCol-1], 18e-3)
    else:
        mle_bounds = (mle_set[iCol-1], mle_set[iCol+1])
    
    # Now create initial guess.
    fileName = mus+f'_cf{cf_set[iRow]:0.1f}Hz_fts{fts:0.2f}_mle{mle_set[iCol]*1e3:0.1f}mm'
    df = df = pd.read_csv(dataDirSim+fileName+'.csv')
    data = df.to_numpy()
    time,lmtc,stim,fsee,_,_ = data.T[0:6]
    t_stimOn, t_stimOff = stimulation.get_stim_timing(time,stim)
    t_stimOn = t_stimOn[-1]-time[-1]
    t_stimOff = t_stimOff[-1]
        
    initialGuess = {}
    initialGuess['cfGuess'] = cf_set[iRow]
    initialGuess['mleGuess'] = mle_set[iCol]
    initialGuess['stimGuess'] = [t_stimOn, t_stimOff]
    
    # Perform optimisation
    AMPO, y, optPar = derive_predictions.opt_ssc_par(2,cf_bounds,fts,mle_bounds,lmtc_avg,muspar,initialGuess)
    time, lmtc, stim, gamma, lcerel, q, lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel = y[0:14]
    
    # Save results
    filename = mus+f'_cfOpt_fts{fts:0.2f}_mleOpt'
    filepath = os.path.join(dataDir,mus,'simsCV',filename+'.csv')
    try:
        df = pd.read_csv(filepath)
        data = df.to_numpy().T
        Pmech = integrate.trapezoid(data[3],data[1])/data[0][-1]
    except:
        Pmech = 0
    
    if AMPO > Pmech:
        # Save data
        data = np.vstack((time,lmtc,stim,fsee,gamma,lcerel,q,fisomrel,fce)).T
        pd.DataFrame(data).to_csv(filepath,index=False,header=['Time [s]','Lmtc [m]',
            'STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]','q [ ]',
            'fisomrel [ ]', 'fce [N]'])

Compute AMPO for imposed MTC length excursion

Code 
"""
This script predicts the maximally attainable AMPO for SSCs with imposed MTC
length excursion using constant MTC shortening/lengthening velocity.

For each imposed MTC length excursion, the script first loads the predictions
for various combinations of cycle frequency and FTS. These predictions are used
to make the initial guess and to define local bounds for the optimisation.
Then, cycle frequency, FTS and stimulation timing are optimised, and the result
is saved to `simsCV`.
"""

#%% Load packages & set directories
import os, sys, pickle
import numpy as np
import pandas as pd
from scipy import integrate
from pathlib import Path

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

import derive_predictions, helpers, stimulation

import matplotlib.pyplot as plt
plt.close('all')

#%% Load muscle parameters
mus = 'GMe3'
parFile = os.path.join(dataDir,mus,mus+'_IM.pkl')
muspar = pickle.load(open(parFile, 'rb'))[0]

# Now we change lsee0 to 3mm
muspar_lsee0_new = 3e-3
eseerelmax = (muspar['fmax']/muspar['ksee'])**0.5/muspar['lsee0']
muspar['ksee'] = muspar['ksee']*(muspar['lsee0']/muspar_lsee0_new)**2
muspar['lsee0'] = muspar_lsee0_new
lmtc_avg = muspar['lce_opt']+(1+eseerelmax)*muspar['lsee0']

#%% Optimise cycle frequency and FTS for imposed MTC length excursion
dataDirSim = os.path.join(dataDir,mus,'simsCV','')

mle_set = np.arange(2, 11.1, 1)*1e-3 # [m]
for mle in mle_set[-1:]:
    # First read out solution for imposed CF, FTS & MLE
    cf_set = np.arange(0.5, 6.1, 0.5) # [Hz]
    fts_set = np.arange(0.05, 0.96, 0.05) # []
    AMPO_set = helpers.load_sims(cf_set,fts_set,mle,mus,dataDirSim)
        
    # Find maximum
    max_idx = np.nanargmax(AMPO_set)
    iRow, iCol = np.unravel_index(max_idx, AMPO_set.shape)
    
    # Create bounds
    try:
        cfBounds = (cf_set[iRow-1], cf_set[iRow+1])
    except:
        cfBounds = (cf_set[iRow-1], 14)
    ftsBounds = (fts_set[iCol-1], fts_set[iCol+1])
    
    # Now create initial guess.
    fileName = mus+f'_cf{cf_set[iRow]:{"0.1f"}}Hz_fts{fts_set[iCol]:{"0.2f"}}_mle{mle*1e3:{"0.1f"}}mm'
    df = df = pd.read_csv(dataDirSim+fileName+'.csv')
    data = df.to_numpy()
    time,lmtc,stim,fsee,_,_ = data.T[0:6]
    t_stimOn, t_stimOff = stimulation.get_stim_timing(time,stim)
    t_stimOn = t_stimOn[-1]-time[-1]
    t_stimOff = t_stimOff[-1]
        
    initialGuess = {}
    initialGuess['cfGuess'] = cf_set[iRow]
    initialGuess['ftsGuess'] = fts_set[iCol]
    initialGuess['stimGuess'] = [t_stimOn, t_stimOff]
    
    # Perform optimisation
    AMPO, y, optPar = derive_predictions.opt_ssc_par(2,cfBounds,ftsBounds,mle,lmtc_avg,muspar,initialGuess)
    time, lmtc, stim, gamma, lcerel, q, lsee, lpee, fisomrel, fsee, fpee, fce, fcerel, vcerel = y[0:14]
    
    # Save results
    filename = mus+f'_cfOpt_ftsOpt_mle{mle*1e3:0.1f}mm'
    filepath = os.path.join(dataDir,mus,'simsCV',filename+'.csv')
    try:
        df = pd.read_csv(filepath)
        data = df.to_numpy().T
        Pmech = integrate.trapezoid(data[3],data[1])/data[0][-1]
    except:
        Pmech = 0
    
    if AMPO > Pmech:
        # Save data
        data = np.vstack((time,lmtc,stim,fsee,gamma,lcerel,q,fisomrel,fce)).T
        pd.DataFrame(data).to_csv(filepath,index=False,header=['Time [s]','Lmtc [m]',
            'STIM [ ]', 'Fsee [N]', 'Gamma [ ]','Lcerel [ ]','q [ ]',
            'fisomrel [ ]', 'fce [N]'])

Analyse simulation results

Code 
"""
This script analyses the constant-velocity SSC predictions after stimulation
timing and SSC parameters have been optimised.

Specifically, the following steps were taken:

-   Load the optimised predictions for all specimens.
-   Interpolate AMPO over a finer cycle frequency, FTS and MTC length excursion
    grid.
-   Average the interpolated AMPO landscape across specimens.
-   Find the parameter combination yielding maximal AMPO.
-   Estimate the 95%-of-maximum ranges for cycle frequency, FTS and MTC length
    excursion.

Custom functions used:

-   `helpers.load_sims(cf_set, fts_set, mle_set, mus, data_dir)`
    :   Load simulation files and compute AMPO values for the requested SSC
        parameter grid.
-   `interpolation.do_4d(data, grid, **kwargs)`
    :   Interpolate the AMPO grid onto a finer parameter grid.
-   `stats.find_max(data, grid)`
    :   Find the maximum value in the interpolated AMPO grid and return the
        corresponding SSC parameter values.
"""

#%% Load packages & set directories
import os, sys, pickle
import numpy as np
from pathlib import Path

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

import helpers, interpolation, stats

#%% Load and interpolate AMPO grids
cfSet   = [0.5, 1, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5, 5.5, 6]
ftsSet  = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95]
mleSet  = [2e-3, 3e-3, 4e-3, 5e-3, 6e-3, 7e-3, 8e-3, 9e-3, 10e-3, 11e-3, 12e-3]

# Better to use a smaller set around optimum
cfSet   = [2.5, 3.0, 3.5, 4.0, 4.5, 5.0]
ftsSet  = [0.70, 0.75, 0.80, 0.85, 0.90, 0.95]
mleSet  = [6e-3, 7e-3, 8e-3, 9e-3, 10e-3, 11e-3]

AMPOsets, AMPOfines, sf = [], [], []
for iMus,mus in enumerate(['GMe1', 'GMe2', 'GMe3']):
    # Load muspar
    parFile = os.path.join(dataDir,mus,mus+'_IM.pkl')
    muspar = pickle.load(open(parFile, 'rb'))[0]
    sf.append(1)
    # sf.append(muspar['lce_opt']*muspar['fmax'])
       
    # Load data   
    dataDir_mus = os.path.join(dataDir,mus,'simsCV')
    AMPOset = helpers.load_sims(cfSet, ftsSet, mleSet,mus,dataDir_mus)
    AMPOfine,(cfFine,ftsFine,mleFine) = interpolation.do_4d(AMPOset,(cfSet,ftsSet,mleSet),N=100,method='linear')    
        
    # Find maximum  
    AMPOset = AMPOset/sf[iMus]
    AMPOfine = AMPOfine/sf[iMus]
    AMPOsets.append(AMPOset)
    AMPOfines.append(AMPOfine)

#%% Average predictions across specimens
AMPOsets = np.stack(AMPOsets, axis=0)
AMPOfines = np.stack(AMPOfines, axis=0)

AMPOmean = np.mean(AMPOfines,0)
AMPOmean = AMPOmean*np.mean(sf)

#%% Find peak AMPO
# Find maximum
AMPOmax, (cfOpt, ftsOpt, mleOpt)  = stats.find_max(AMPOmean,(cfFine,ftsFine, mleFine))

print("AMPO = %1.2f mW" % (AMPOmax*1e3))
print("CF = %1.2f Hz" % cfOpt)
print("FTS = %1.2f" % ftsOpt)
print("MLE = %1.2f mm" % (mleOpt*1e3))

#%% Estimate 95%-of-maximum ranges
iRow, iCol, iDep = np.unravel_index(np.nanargmax(AMPOmean), AMPOmean.shape)

maxAMPO = np.nanmax(AMPOmean)

# 95%
AMPO_CF = AMPOmean[:,iRow,iDep] # @opt cf
idx = np.where(AMPO_CF>0.95*maxAMPO)
cf95 = cfFine[idx]
cfRange = (cf95[0], cf95[-1])
print(cfRange)

AMPO_FTS = AMPOmean[iRow,:,iDep] # @opt fts
idx = np.where(AMPO_FTS>0.95*maxAMPO)
fts95 = ftsFine[idx]
ftsRange = (fts95[0], fts95[-1])
print(ftsRange)

AMPO_AMP = AMPOmean[iRow,iCol,:] # @opt amp
idx = np.where(AMPO_AMP>0.95*maxAMPO)
mle95 = mleFine[idx]
mleRange = (mle95[0]*1e3, mle95[-1]*1e3)
print(mleRange)