<jupyter_start><jupyter_text>Empirical method to estimate ke and ki from experimental stationary distributions<jupyter_code>#from sys import platform
import os
# rSNAPsim Image Processing
cwd = os.getcwd() # current working directory
os.chdir('..') # move one directory up
os.chdir('..') # move one directory up
import rSNAPsim_IP as rss_IP #importing rSNAPsim.
os.chdir(cwd) # return to the working directory<jupyter_output>_____ _   _____    ____       _              ________ 
   _____/ ___// | / /   |  / __ \_____(_)___ ___     /  _/ __ \
  / ___/\__ \/  |/ / /| | / /_/ / ___/ / __ `__ \    / // /_/ /
 / /   ___/ / /|  / ___ |/ ____(__  ) / / / / / /  _/ // ____/ 
/_/   /____/_/ |_/_/  |_/_/   /____/_/_/ /_/ /_/  /___/_/<jupyter_text>Simulation Workflow   Is it possible to go in the  inverse   direction?   Testing if each combination of {ki,ke,L}  maps  to a specific distribution  {$\mu$,$\sigma$}.   The goal is to find  phenomenological equations  $f_{ke}(\mu,\sigma,L)$ and $f_{ki}(\mu,\sigma,L)$<jupyter_code># Importing libraries

from skimage import io ; from skimage.io import imread; from skimage.measure import find_contours
# Plotting
import matplotlib as mpl ; import matplotlib.pyplot as plt 
# To manipulate arrays
import numpy as np 
# To create interactive elements
#import ipywidgets as widgets ; from ipywidgets import interact, interactive, fixed, interact_manual, Button, HBox, VBox, Layout, GridspecLayout ; from IPython.display import Image
# Figure style and size
mpl.rc('figure',  figsize=(7, 7)); mpl.rc('image', cmap='gray'); mpl.rc('axes', grid=False)
# change the following to %matplotlib notebook for interactive plotting
import tifffile
import os; from os import listdir; from os.path import isfile, join
import random
# Parallel computing
from joblib import Parallel, delayed
import multiprocessing
from tqdm import tqdm
# Optimization
import scipy as sp
from scipy.optimize import curve_fit
# 3D plotting
from mpl_toolkits import mplot3d

# Deffining directory
save_to_path = '/home/luisub/Desktop/image_processing_toolbox/Output/testing_distributions'
rSNAPsim_directory = '/home/luisub/Documents/rSNAPsim-dev' 
if not os.path.exists(save_to_path):
    os.makedirs(save_to_path)
##################################
##################################
## Stochastic Simulations
number_spots = 100
n_frames = 200 # in seconds
step_size = 1 # in seconds
starting_time = 1000
##################################
##################################
#### Section that runs the stochastic simulations
sequences_directory = rSNAPsim_directory+'/gene_files/Bactin_withTags.txt'
## RUNNING THE SSA
# current working directory
cwd = os.getcwd()
#importing rSNAPsim.
os.chdir(rSNAPsim_directory)    
import rSNAPsim
# rename rSNAPsim as rss
rss = rSNAPsim.rSNAPsim()
# returning to the directory with the image processing codes.
os.chdir(cwd)    
rss.open_seq_file(sequences_directory)
# now to look for open reading frames
rss.get_orfs(rss.sequence_str)
rss.orfs# reading the 3 ORF.
# Deffining a Protein passing the DNA sequence 
protein1 = rss.sequence_str  
rss.nt2aa(protein1) # Converting nucleotides to Protein
rss.get_temporal_proteins()  #proteins sorted by
# Analyze the first protein, analyze function takes the aa_seq then nt_seq
rss.analyze_poi(rss.pois[0],rss.pois_seq[0]) 
rss.POI
rss.POI.total_length
rss.POI.tag_epitopes<jupyter_output><empty_output><jupyter_text>SSA for $\beta$-actin witn different ke and ki  Simulaitng multiple ke and (ki=0.03)<jupyter_code>number_of_evaluations_short = 5
list_simulations_ke = []
ki = 0.03
list_ke_short = np.linspace(1,number_of_evaluations_short,number_of_evaluations_short)
for i in range (0,len(list_ke_short)):
    ssa = rss.ssa_solver(n_traj = number_spots, start_time=starting_time,tf=starting_time+n_frames, tstep=starting_time+n_frames,k_elong_mean=list_ke_short[i], k_initiation=ki)  # tstep = total number of steps including the burnin time 
    simulated_trajectories_ke = ssa.intensity_vec
    list_simulations_ke.append (simulated_trajectories_ke)
# Plotting
fig, axes = plt.subplots(nrows=1, ncols=len(list_ke_short), figsize=(15, 5))
for i, ke_val in enumerate (list_ke_short):
    axes[i].hist(list_simulations_ke[i][:,-1], alpha=0.5,bins=30,density=True, stacked=True, color=(random.random(),random.random(),random.random()) )      
    axes[i].set(title= 'ke = ' + str(round(ke_val,1)))
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Simulaitng multiple ki and (ke = 3)<jupyter_code>list_simulations_ki = []
ke = 3
list_ki_short = np.linspace(0.02,0.08,number_of_evaluations_short)
for i in range (0,len(list_ki_short)):
    ssa = rss.ssa_solver(n_traj = number_spots, start_time=starting_time,tf=starting_time+n_frames, tstep=starting_time+n_frames,k_elong_mean=ke, k_initiation=list_ki_short[i])  # tstep = total number of steps including the burnin time 
    simulated_trajectories_ki = ssa.intensity_vec
    list_simulations_ki.append (simulated_trajectories_ki)
# Plotting
fig, axes = plt.subplots(nrows=1, ncols=len(list_simulations_ki), figsize=(15, 5))
for i, ki_val in enumerate (list_ki_short):
    axes[i].hist(list_simulations_ki[i][:,-1], alpha=0.5,bins=30,density=True, stacked=True, color=(random.random(),random.random(),random.random()) )      
    axes[i].set(title= 'ki = ' + str(round(ki_val,2)))
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Changing ke and ki simultaneously<jupyter_code>number_of_evaluations = 10
list_ke = np.linspace(1,10,number_of_evaluations)     # np.linspace(min, max, n)
list_ki = np.linspace(0.02,0.08,number_of_evaluations)
def run_index (ki,ke):
    return ke,ki#, std_ssa
matrix_index = [[run_index(list_ki[j],list_ke[i]) for i in range(number_of_evaluations)] for j in range(number_of_evaluations)] 
ind_array = np.asarray(matrix_index)
#print('ke values: ')
#print(ind_array[:,:,0])
#print('ki values: ')
#print(np.round(ind_array[:,:,1],2))
list_simulations_ke[0].shape
def run_ssa (ki,ke):
    try:
        ssa = rss.ssa_solver(n_traj = number_spots, start_time=starting_time,tf=starting_time+n_frames, tstep=starting_time+n_frames,k_elong_mean=ke, k_initiation=ki)  # tstep = total number of steps including the burnin time 
        simulated_trajectories = ssa.intensity_vec
        mean_ssa = np.mean(simulated_trajectories[:,-1])
        std_ssa = np.std(simulated_trajectories[:,-1])
    except:
        mean_ssa = 0
        std_ssa = 0
    return round(mean_ssa,2), round(std_ssa,2)<jupyter_output><empty_output><jupyter_text>Loop that iterates for multiple ki and ke<jupyter_code>ssa_matrix = [[run_ssa(list_ki[j],list_ke[i]) for i in range(number_of_evaluations)] for j in range(number_of_evaluations)]
ssa_values=np.array(ssa_matrix) # converting the list into numpy array
means =np.round(ssa_values[:,:,0],0) # mean values
stds =np.round(ssa_values[:,:,1],0) # standard deviations
means
stds
# Plotting means
fig, ax = plt.subplots( figsize=(5, 5))
pos = ax.imshow(means,origin='lower', cmap='Reds', interpolation='none')
ax.set(title= 'Mean values')
ax.set_xlabel('ke')
ax.set_ylabel('ki')
cbar = fig.colorbar(pos, ax=ax)
cbar.ax.set_title('Mean \n', fontsize=9)
plt.show()
# Plotting stdvs
fig, ax = plt.subplots( figsize=(5, 5))
pos = ax.imshow(stds,origin='lower', cmap='Blues', interpolation='none')
ax.set(title= 'Std values')
ax.set_xlabel('ke')
ax.set_ylabel('ki')
cbar = fig.colorbar(pos, ax=ax)
cbar.ax.set_title('Std \n', fontsize=9)
plt.show()<jupyter_output><empty_output><jupyter_text>It appears that each combination of {ki,ke,L}  maps  to a specific distribution {$\mu$,$\sigma$}.<jupyter_code>X, Y = np.meshgrid(list_ke, list_ki )
fig = plt.figure(figsize=(10, 5))
ax = plt.axes(projection='3d')
# fourth dimention - colormap
color_dimension = stds # change to desired fourth dimension
minn, maxx = color_dimension.min(), color_dimension.max()
norm = mpl.colors.Normalize(minn, maxx)
m = plt.cm.ScalarMappable(norm=norm, cmap='viridis')
m.set_array([])
fcolors = m.to_rgba(color_dimension)
surf =ax.plot_surface(X, Y, means, rstride=1, cstride=1, facecolors=fcolors, vmin=minn, vmax=maxx, shade=False)
#ax.view_init(60, 35)
ax.set_xlabel('ke')
ax.set_ylabel('ki')
ax.set_zlabel('mean')
cbar = fig.colorbar(m, shrink=0.6, aspect=10)
cbar.ax.set_title('Std \n', fontsize=9)
plt.show()<jupyter_output><empty_output><jupyter_text>Testing if the distributions follow a pattern<jupyter_code>def exp_decay(x, a, b, c):
     return a*np.exp(-b*x) + c<jupyter_output><empty_output><jupyter_text>Mean values  \begin{equation*}\mu(ke)   = A \cdot e^{-B*ke} + C \end{equation*}  Mean values (fixed ki = 0.02)<jupyter_code># Optimization for the means
popt_mu, pcov_mu = curve_fit(exp_decay, list_ke, means[0,:]) # optimization
fit_mu = exp_decay(list_ke, popt_mu[0],popt_mu[1],popt_mu[2]) # generating line with optimized values

# Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
ax.plot(list_ke,means[0,:], 'ok',linewidth = 3,label='Data') # mean
ax.plot(list_ke,fit_mu, '-r', markersize=10,label=' $f(ke) = %.0f * e^{-%.1f*ke} + %.1f$' % (popt_mu[0],popt_mu[1],popt_mu[2])) # mean
ax.set(title= 'Mean values')
ax.set_xlabel('ke')
ax.set_ylabel('Mean')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Showing for all parameters<jupyter_code># Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
for i in range(number_of_evaluations):
    ax.plot( list_ke, means[i,:] ,label='ki = ' + str(round(list_ki[i],2) )) # mean
ax.set(title= 'Mean values')
ax.set_xlabel('ke')
ax.set_ylabel('Mean')
ax.legend()
plt.show()
# Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
for i in range(number_of_evaluations):
    ax.plot( list_ke, means[i,:]/means[i,0],label='ki = ' + str(round(list_ki[i],2) )) # mean
ax.set(title= 'Normalized Mean values')
ax.set_xlabel('ke')
ax.set_ylabel('Mean/Mean[0]')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Linear function for the inital mean value \begin{equation*}A = \mu[0]   = m\cdot ki + b\end{equation*}<jupyter_code>def linear(x, m, b):
     return m*x + b
popt_init_mean, pcov_init_mean = curve_fit(linear, list_ki, means[:,0])
fit_init_mean= linear(list_ki, popt_init_mean[0],popt_init_mean[1])

# Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
ax.plot(list_ki, means[:,0], 'ok',linewidth = 3,label='Data') # mean
ax.plot(list_ki,fit_init_mean, '-r', markersize=10,label='$f(ki) = %.1f*x + %.1f$' % (popt_init_mean[0],popt_init_mean[1])) # mean
ax.set(title= 'A = Mean [0]')
ax.set_xlabel('ki')
ax.set_ylabel('A = Mean [0]')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>A phenomenological function for $f_{\mu}(ke,ki)$ \begin{equation*}\mu(ke,ki)   = ki (c_{1} \cdot e^{-B*ke} + c_{2})\end{equation*}   Standard Deviation   Standard deviation values (fixed ki = 0.02)<jupyter_code># Optimization for the standar deviation
popt_sd, pcov_mu = curve_fit(exp_decay, list_ke, stds[0,:]) # optimization
fit_sd = exp_decay(list_ke, popt_sd[0],popt_sd[1],popt_sd[2]) # generating line with optimized values

# Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
ax.plot(list_ke,stds[0,:], 'ok',linewidth = 3,label='Data') # mean
ax.plot(list_ke,fit_sd, '-r', markersize=10,label='$f(mu) = %.0f * e^{-%.1f*ke} + %.1f$' % (popt_sd[0],popt_sd[1],popt_sd[2])) # mean
ax.set(title= 'Standard deviation values')
ax.set_xlabel('ke')
ax.set_ylabel('Standard deviation')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Showing for all parameters<jupyter_code># Plotting the standar deviation
fig, ax = plt.subplots( figsize=(10, 5))
for i in range(number_of_evaluations):
    ax.plot( list_ke, stds[i,:]  ,label='ki = ' + str(round(list_ki[i],2) )) # mean
ax.set(title= 'Std values')
ax.set_xlabel('ke')
ax.set_ylabel('std')
ax.legend()
plt.show()
# Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
for i in range(number_of_evaluations):
    ax.plot( list_ke, stds[i,:]-np.amin(stds[i,:]),label='ki = ' + str(round(list_ki[i],2) )) # mean
ax.set(title= 'Normalized std values')
ax.set_xlabel('ke')
ax.set_ylabel('std-std_min')
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>A phenomenological function for $f_{\sigma}(ke,ki)$ \begin{equation*}\sigma(ke,ki)   = ki (\hat{c_{1}} \cdot e^{-\hat{B}*ke} + \hat{c_{2}})\end{equation*}  The goal is to find  phenomenological equations  $f_{ke}(\mu,\sigma,L)$ and $f_{ki}(\mu,\sigma,L)$  Standard deviation as a function of mean<jupyter_code>def hill(x, vmax, km,n):
     return vmax*(x**n/(km+x**n))<jupyter_output><empty_output><jupyter_text>Hill function \begin{equation*}y   = vmax \frac{x^{n}} {km + x^{n}}\end{equation*}<jupyter_code>p0 = np.array([30,60,1])
popt_combined, pcov_combined = curve_fit(hill, fit_mu, fit_sd,p0)
fit_combined= hill(fit_mu, popt_combined[0],popt_combined[1],popt_combined[2])

# Plotting the means
fig, ax = plt.subplots( figsize=(10, 5))
for i in range(0,number_of_evaluations):
    ax.plot(means[i,:], ((stds[i,:]**2)/10)/factor, '-',linewidth = 3,label='ki = ' + str(round(list_ki[i],2) )) 

plt.plot([0,100],[0,100],'-k')    
#ax.plot(fit_mu,fit_combined, '-r', markersize=10,label='$f(mu) = %.0f * (mu^{%.2f} / %.0f + mu^{%.2f} )$' % (popt_combined[0],popt_combined[2],popt_combined[1],popt_combined[2])) # mean
ax.set(title= 'Mean_vs_Std')
ax.set_xlabel('Mean (mu)')
ax.set_ylabel('Standard deviation')
ax.legend()
plt.show()
stds[0,:]
means[0,:]
np.sqrt(stds[0,:])<jupyter_output><empty_output><jupyter_text>Obtaning ke from experimental $\mu$\begin{equation*}ke(\mu)   = \frac{ln(A) -ln(\mu - C)} {B}\end{equation*}<jupyter_code>factor  = (1-(318/(2*712))) / (1-((2*318/(3*712))) ) 
factor<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16705641625255774, 'fraction_numerical': 0.0176118094124525, 'mean_word_length': 5.939655172413793, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '477065', 'license': 'MIT', 'n_tokens_mistral': 5503, 'n_tokens_neox': 5272, 'n_words': 1031, 'path': 'rsnapsim/interactive_notebooks/Empirical_estimation_mu_sigma.ipynb', 'repo_name': 'MunskyGroup/rSNAPsim', 'size': '13684'}