Simulated data

On this page you can find the code to make the simulated data described in the section: “Simulated data”.

Obtain parameters

The following simple script was used to create pickles with the ‘actual’ parameter values, which were obtained from literature (see manuscript).

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

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

# %% Make dict with parameters values and save as pickle
for mus in ['GMs1','GMs2','GMs3']:
    muspar = {}
    
    # Force-length relationship
    muspar['w'] = 0.5;
    muspar['n'] = 2
    
    # Force-velocity relationship
    muspar['fasymp'] = 1.5
    muspar['slopfac'] = 2
    muspar['vfactmin'] = 0.1
    
    # Activation dynamics    
    muspar['a_act'] = -7.369163055099003
    muspar['b_act'] = np.array([5.170927028993413, 0.5955111970420514, 0])
    muspar['q0'] = 0.005
    muspar['gamma_0'] = 1e-5
    muspar['kCa'] = 0.8e-5
    muspar['tact'] = 1/37
    muspar['tdeact'] = 1/37
      
    fileName = os.path.join(dataDir,'GM_muspar.pkl')
    pickle.dump(muspar, open(fileName, 'wb'))
    
    # Rat specific
    if mus == 'GMs1':
        muspar['lce_opt'] = 1.32e-2 # [m]
        muspar['fmax'] = 13.39 # [N]
        muspar['ksee'] = 4.22e6
        muspar['kpee'] = 2.13e5
        muspar['lsee0'] = 2.83e-2 # [m]
        muspar['lpee0'] = 1.39e-2 # [m]
        # muspar['eseerelmax'] = (muspar['fmax']/4.22e6)**(1/2)/muspar['lsee0']
        # muspar['epeerelmax'] = (muspar['fmax']/2.13e5)**(1/2)/muspar['lpee0']
        muspar['a'] = 0.20*muspar['fmax']
        muspar['b'] = 3.15*muspar['lce_opt']
    elif mus == 'GMs2':
        muspar['lce_opt'] = 1.23e-2 # [m]
        muspar['fmax'] = 13.81 # [N]
        muspar['ksee'] = 3.64e6
        muspar['kpee'] = 1.65e5
        muspar['lsee0'] = 3.03e-2 # [m]
        muspar['lpee0'] = 1.32e-2 # [m]
        # muspar['eseerelmax'] = (muspar['fmax']/3.64e6)**(1/2)/muspar['lsee0']
        # muspar['epeerelmax'] = (muspar['fmax']/1.65e5)**(1/2)/muspar['lpee0']
        muspar['a'] = 0.13*muspar['fmax']
        muspar['b'] = 2.02*muspar['lce_opt']
    elif mus == 'GMs3':
        muspar['lce_opt'] = 1.12e-2 # [m]
        muspar['fmax'] = 12.28 # [N]
        muspar['ksee'] = 3.47e6
        muspar['kpee'] = 5.11e5
        muspar['lsee0'] = 2.65e-2 # [m]
        muspar['lpee0'] = 1.34e-2 # [m]
        # muspar['eseerelmax'] = (muspar['fmax']/3.47e6)**(1/2)/muspar['lsee0']
        # muspar['epeerelmax'] = (muspar['fmax']/5.11e5)**(1/2)/muspar['lpee0']
        muspar['a'] = 0.21*muspar['fmax']
        muspar['b'] = 3.73*muspar['lce_opt']
    else:
        print('Incorrect muscle selected!')
    
    filepath = os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl')
    pickle.dump(muspar, open(filepath, 'wb'))

Quick-release data

#%% Load modules
import pickle, os, sys, glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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 make_data, stats

plt.close('all')

#%% Readout muscle parameter values
if not 'mus' in locals():
    mus = 'GMs1'

parFile = os.path.join(dataDir,mus,'parameters',mus+'_OR.pkl')
muspar = pickle.load(open(parFile, 'rb'))
eseeMax = (muspar['fmax']/muspar['ksee'])**(1/2) # [mm] SEE elongation @ fmax
lmtcOpt = eseeMax + muspar['lsee0']+muspar['lce_opt']

#%% Simulate QR data - 'Experimental data'
# Setting which are similar for each QR
optsQR = {
    'durIso':   [0.3, 0.2], # [s] duration of 1st and 2nd isometric phase
    'durStep':  10e-3,      # [s] duration of step
    'dLmtc':    -0.2e-3,     # [m] MTC length change during step 
    'tStim':    [0.1, 0.4]  # [s] stim onset and offset time
    }
dataDirQR = os.path.join(dataDir,mus,'dataExp','QR','')

# Define over which range we will simulate QR
lmtcRange = np.linspace(stats.floor(lmtcOpt,3)+3e-3,stats.floor(lmtcOpt,3)-6e-3,10)

# First delete all files in folder
files = glob.glob(dataDirQR+r'\*')    
for iFile,filename in enumerate(files):
    os.remove(filename)

# Make & save data
plt.close('all')
fig, ax = plt.subplots(3,1)
for iFile,lmtc0 in enumerate(lmtcRange):
    data = make_data.qr(lmtc0,muspar,optsQR)
    
    # Check
    ax[0].plot(data[:,0],data[:,2]) # sitm
    ax[1].plot(data[:,0],data[:,1]) # lmtc
    ax[2].plot(data[:,0],data[:,3]) # fsee
    
    fileName = mus+'_QR'+'{:02d}'.format(iFile+1)+'_OR.csv'
    pd.DataFrame(data).to_csv(dataDirQR+fileName,index=False, 
                                header=['time [s]','Lmtc [m]','STIM [ ]','Fsee [N]'])
    print('QR: '+mus+', '+'OR'+', File: '+str(iFile+1))

#%% Simulate QR data - Monte Carlo simulations
for iPar in range(1,51):
    dataDirQR = os.path.join(dataDir,mus,'dataMC','QR',f'{iPar:{"02d"}}','')
    
    # Define over which range we will simulate QR, this is a different range (i.e., slightly shifted). 
    # This is what one will (or should do) in experiment too!
    np.random.seed(iPar)
    shifts = np.random.rand(10)*1e-3
    
    # First delete all files in folder
    files = glob.glob(dataDirQR+r'\*')    
    for iFile,filename in enumerate(files):
        os.remove(filename)
    
    # Make & save data
    plt.close('all')
    fig, ax = plt.subplots(3,1)
    for iFile,lmtc0 in enumerate(lmtcRange):
        mcpar = muspar.copy()
        # Distrub data
        eseeMax = lmtcOpt+shifts[iFile]-muspar['lce_opt']-muspar['lsee0'] # [mm] SEE elongation @ fmax
        mcpar['ksee'] = mcpar['fmax']/(eseeMax)**2
        data = make_data.qr(lmtc0,mcpar,optsQR)
        
        # Check
        ax[0].plot(data[:,0],data[:,2]) # sitm
        ax[1].plot(data[:,0],data[:,1]*1e3) # lmtc
        ax[2].plot(data[:,0],data[:,3]) # fsee
        
        fileName = mus+'_QR'+'{:02d}'.format(iFile+1)+'_MC{:02d}'.format(iPar)+'.csv'
        pd.DataFrame(data).to_csv(dataDirQR+fileName,index=False, 
                                  header=['time [s]','Lmtc [m]','STIM [ ]','Fsee [N]'])
        print('QR: '+mus+', '+'MC{:02d}'.format(iPar)+', File: '+str(iFile+1))

Step-ramp data

#%% Load modules
import pickle, os, sys, glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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 make_data, hillmodel, stats

plt.close('all')

#%% Readout muscle parameter values
if not 'mus' in locals():
    mus = 'GMs1'

parFile = os.path.join(dataDir,mus,'Parameters',mus+'_OR.pkl')
muspar = pickle.load(open(parFile, 'rb'))
eseeMax = (muspar['fmax']/muspar['ksee'])**(1/2) # [mm] SEE elongation @ fmax
lmtcOpt = eseeMax + muspar['lsee0']+muspar['lce_opt']
lmtc0 = lmtcOpt+0.5e-3 # [m] MTC-length at t=0

#%% Simulate SR data - 'Experimental data'
fceRange = np.array([0.85, 0.72, 0.60, 0.49, 0.38, 0.28, 0.20, 0.14, 0.09])*muspar['fmax']

# Setting which are similar for each SR
optsSR = {
    'durIso':   0.3,        # [s] duration of 1st isometric phase
    'durStep':  10e-3,      # [s] duration of step
    'durRamp':  0.2,        # [m] MTC length change during step
    'tStim':    [0.1, 0.4]  # [s] stim onset and offset time
    }
dataDirSR = os.path.join(dataDir,mus,'dataExp','SR','')

# First delete all files in folder
files = glob.glob(dataDirSR+r'\*')    
for iFile,filename in enumerate(files):
    os.remove(filename)

# Make & save data
plt.close('all')
fig, ax = plt.subplots(3,1)
for iFile,fce in enumerate(fceRange):
    vRamp = hillmodel.Fce2Vce(fce,1,1,muspar)[0]
    esee = (fceRange[iFile]/muspar['ksee'])**0.5
    eseeMax = (muspar['fmax']/muspar['ksee'])**0.5
    dLmtc = np.round((esee-eseeMax)+(vRamp*1e-2)*0.5,5)
    data = make_data.sr(lmtc0,dLmtc,vRamp,muspar,optsSR)
    
    # Check
    ax[0].plot(data[:,0],data[:,2]) # sitm
    ax[1].plot(data[:,0],data[:,1]) # lmtc
    ax[2].plot(data[:,0],data[:,3]) # fsee
    fileName = mus+'_SR'+'{:02d}'.format(iFile+1)+'_OR.csv'
    pd.DataFrame(data).to_csv(dataDirSR+fileName,index=False, 
                              header=['time [s]','Lmtc [m]','STIM [ ]','Fsee [N]'])
    print('SR: '+mus+', '+'OR'+', File: '+str(iFile+1))

#%% Simulate SR data - Monte Carlo simulations
for iPar in range(1,51):
    dataDirSR = os.path.join(dataDir,mus,'dataMC',f'{iPar:{"02d"}}','SR','')
    
    # Define MTC shifts, we will not choose a different LMTC length at the start, 
    # because it's unrealistic that one will change this for every SR in an experiment
    np.random.seed(iPar)
    shifts = np.random.rand(9)*1e-3
    
    # First delete all files in folder
    files = glob.glob(dataDirSR+r'\*')    
    for iFile,filename in enumerate(files):
        os.remove(filename)
    
    # Make & save data
    plt.close('all')
    fig, ax = plt.subplots(3,1)
    for iFile,fce in enumerate(fceRange):
        mcpar = muspar.copy()
        # Distrub data
        eseeMax = lmtcOpt+shifts[iFile]-muspar['lce_opt']-muspar['lsee0'] # [mm] SEE elongation @ fmax
        mcpar['ksee'] = mcpar['fmax']/(eseeMax)**2
    
        vRamp = hillmodel.Fce2Vce(fceRange[iFile],1,1,mcpar)[0]
        esee = (fceRange[iFile]/mcpar['ksee'])**0.5
        eseeMax = (mcpar['fmax']/mcpar['ksee'])**0.5
        dLmtc = np.round((esee-eseeMax)+(vRamp*1e-2)*0.5,5)
        data = make_data.sr(lmtc0,dLmtc,vRamp,mcpar,optsSR)
        
        # Check
        ax[0].plot(data[:,0],data[:,2]) # sitm
        ax[1].plot(data[:,0],data[:,1]*1e3) # lmtc
        ax[2].plot(data[:,0],data[:,3]) # fsee
        
        fileName = mus+'_SR'+'{:02d}'.format(iFile+1)+'_MC{:02d}'.format(iPar)+'.csv'
        pd.DataFrame(data).to_csv(dataDirSR+fileName,index=False, 
                                  header=['time [s]','Lmtc [m]','STIM [ ]','Fsee [N]'])
        print('SR: '+mus+', '+'MC{:02d}'.format(iPar)+', File: '+str(iFile+1))

Isometric data

#%% Load modules
import pickle, os, sys, glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
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 make_data, stats

plt.close('all')

#%% Readout muscle parameter values
if not 'mus' in locals():
    mus = 'GMs3'

parFile = os.path.join(dataDir,mus,'Parameters',mus+'_OR.pkl')
muspar = pickle.load(open(parFile, 'rb'))
eseeMax = (muspar['fmax']/muspar['ksee'])**(1/2) # [mm] SEE elongation @ fmax
lmtcOpt = eseeMax + muspar['lsee0']+muspar['lce_opt']

#%% Generate ACT data
if mus == 'GMs1':
    lmtcRange = np.array([43,41,39,37, 43,41,39,37, 43,41,39,37])*1e-3 # [m] isometric MTC length
elif mus == 'GMs2':
    lmtcRange = np.array([44,42,40,38, 44,42,40,38, 44,42,40,38])*1e-3 # [m] isometric MTC length
elif mus == 'GMs3':
    lmtcRange = np.array([39,37,35,33, 39,37,35,33, 39,37,35,33])*1e-3 # [m] isometric MTC length
# durStimRange = np.array([5,5,5, 10,10,10, 15,15,15])*1e-3 # [s] stimulation duration
durStimRange = np.array([35,35,35,35 , 65,65,65, 65, 95,95,95,95])*1e-3 # [s] stimulation duration

#%% Simulate ISOM data - 'Experimental data'
dataDirISOM = os.path.join(dataDir,mus,'dataExp','ISOM','')

# First delete all files in folder
files = glob.glob(dataDirISOM+r'\*')    
for iFile,filename in enumerate(files):
    os.remove(filename)

# Make & save data
plt.close('all')
fig, ax = plt.subplots(3,4)
for iFile,(lmtc0,durStim) in enumerate(zip(lmtcRange,durStimRange)):
    data = make_data.isom(lmtc0,durStim,muspar)
    
    # Check
    iT = np.mod(iFile,4)
    ax[0,iT].plot(data[:,0],data[:,2]) # sitm
    ax[1,iT].plot(data[:,0],data[:,1]) # lmtc
    ax[2,iT].plot(data[:,0],data[:,3]) # fsee
    
    fileName = mus+'_ISOM'+'{:02d}'.format(iFile+1)+'_OR.csv'
    pd.DataFrame(data).to_csv(dataDirISOM+fileName,index=False, 
                                header=['time [s]','Lmtc [m]','STIM [ ]','Fsee [N]'])
    print('ISOM: '+mus+', '+'OR'+', File: '+str(iFile+1))

#%% Simulate ISOM data - Monte Carlo simulations
for iPar in range(1,51):
    dataDirISOM = os.path.join(dataDir,mus,'dataMC',f'{iPar:{"02d"}}','ISOM','')
    
    # Define MTC shifts, we will not choose a different LMTC length at the start, 
    # because it's unrealistic that one will change this for every SR in an experiment
    np.random.seed(iPar)
    shifts = np.random.rand(12)*1e-3
    
    # First delete all files in folder
    files = glob.glob(dataDirISOM+r'\*')    
    for iFile,filename in enumerate(files):
        os.remove(filename)
    
    # Make & save data
    plt.close('all')
    fig, ax = plt.subplots(3,4)
    for iFile,(lmtc0,durStim) in enumerate(zip(lmtcRange,durStimRange)):
        mcpar = muspar.copy()
        # Distrub data
        eseeMax = lmtcOpt+shifts[iFile]-muspar['lce_opt']-muspar['lsee0'] # [mm] SEE elongation @ fmax
        mcpar['ksee'] = mcpar['fmax']/(eseeMax)**2
        data = make_data.isom(lmtc0,durStim,mcpar)
        
        # Check
        iT = np.mod(iFile,4)
        ax[0,iT].plot(data[:,0],data[:,2]) # sitm
        ax[1,iT].plot(data[:,0],data[:,1]) # lmtc
        ax[2,iT].plot(data[:,0],data[:,3]) # fsee
        
        fileName = mus+'_ISOM'+'{:02d}'.format(iFile+1)+'_MC{:02d}'.format(iPar)+'.csv'
        pd.DataFrame(data).to_csv(dataDirISOM+fileName,index=False, 
                                  header=['time [s]','Lmtc [m]','STIM [ ]','Fsee [N]'])
        print('ISOM: '+mus+', '+'MC{:02d}'.format(iPar)+', File: '+str(iFile+1))