fig1_functions_latest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fig 1 just functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import my packages\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "from mpl_toolkits import mplot3d\n",
    "from matplotlib import animation\n",
    "from scipy.fft import fft, ifft, ifft2, fftfreq, fft2, fftshift\n",
    "from scipy import signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_FT(FWHM, x_range): # y_range = 0 1D, y_range = x_range 2D\n",
    "    x = x_range\n",
    "    sigma = FWHM*(1/(2*(2*np.log(2))**(1/2))) \n",
    "    h = (1 / (2 * np.pi * (sigma ** 2)) ** (1 / 2)) * np.exp(-(x ** 2) / (2 * sigma ** 2))\n",
    "    a = fft(h)\n",
    "    sample_freq_x = fftfreq(len(x), d=0.5)\n",
    "    h_shifted = fftshift(a)\n",
    "    sample_freq_x = fftshift(sample_freq_x)\n",
    "    return(h,h_shifted,sample_freq_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### temporal filters and their FTs\n",
    "def time_filter(tau,t, time_step):\n",
    "    theta = 1 # Define Heaviside for t-t'>0 , t'=0\n",
    "    f = 2*tau**(-3/2)*t*theta*np.exp(-t/tau)\n",
    "    f_dev = 2*tau**(-3/2)*(tau-t)*theta*np.exp(-t/tau)\n",
    "    a=fft(f)\n",
    "    amplitude = np.abs(a)\n",
    "    n= len(t)\n",
    "    sample_freq = fftfreq(n, time_step)\n",
    "    # Shift the Fourier transform to center the zero-frequency component\n",
    "    f_shifted = fftshift(a)\n",
    "    sample_freq_shifted = fftshift(sample_freq)\n",
    "    return(sample_freq_shifted,f_shifted)\n",
    "def time_filter_der(tau, t, time_step):\n",
    "    theta = 1 # Define Heaviside for t-t'>0 , t'=0\n",
    "    f = 2*tau**(-3/2)*t*theta*np.exp(-t/tau)\n",
    "    f_dev = 2*tau**(-3/2)*(tau-t)*theta*np.exp(-t/tau)\n",
    "    a=fft(f_dev)\n",
    "    amplitude = np.abs(a)\n",
    "    n= len(t)\n",
    "    sample_freq = fftfreq(n, time_step)\n",
    "    # Shift the Fourier transform to center the zero-frequency component\n",
    "    f_dev_shifted = fftshift(a)\n",
    "    sample_freq_shifted_dev = fftshift(sample_freq)\n",
    "    return(sample_freq_shifted_dev,f_dev_shifted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def edge_stimulus(which_stimulus, t_range, x_range, x_initial): #FTB_ON,FTB_OFF,BTF_OFF,BTF_ON #t = np.arange(0, 8, 1/240) # x = np.arange(-180, 180, 0.5)\n",
    "    if which_stimulus == 'FTB_ON':\n",
    "        dx = -(x_range[:, np.newaxis] - 30 * t_range - x_initial)\n",
    "        c = np.heaviside(dx, 0)\n",
    "    elif which_stimulus == 'FTB_OFF':\n",
    "        dx = (x_range[:, np.newaxis] - 30 * t_range - x_initial)\n",
    "        c = np.heaviside(dx, 0)\n",
    "    elif which_stimulus == 'BTF_OFF':\n",
    "        dx = (-x_range[:, np.newaxis] - 30 * t_range - x_initial)\n",
    "        c = np.heaviside(dx, 0)\n",
    "    elif which_stimulus == 'BTF_ON':\n",
    "        dx = -(-x_range[:, np.newaxis] - 30 * t_range - x_initial)\n",
    "        c = np.heaviside(dx, 0)\n",
    "    return(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gratings(c_o, f, k, t_range, x_range):\n",
    "    omega=2*math.pi*f\n",
    "    t=t_range \n",
    "    x=x_range # x in degrees, I just need x to be larger than  t\n",
    "    contrast = np.empty([len(x),len(t)]) # now instead of creating empty arrays with contrast = np.array([]) and fill with append, I declare the size & it's faster\n",
    "    for j, x_value in enumerate(x):\n",
    "        contrast[j,:] = c_o*np.sin(omega*t-k*x_value)   # c= c_o*np.sin(omega*t[i]-k*x[k]) \n",
    "    c = contrast\n",
    "    return(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contrast_fourier(c,x_range,t_range):\n",
    "    # contrast signal Fourier space\n",
    "    t = t_range\n",
    "    x = x_range\n",
    "    a=fft2(c)\n",
    "    c_shifted = fftshift(a)\n",
    "    #amplitude = np.abs(a)\n",
    "    # Calculate the corresponding spatial and temporal frequencies\n",
    "    spatial_freq = fftfreq(len(x), d=x[1] - x[0])\n",
    "    temporal_freq = fftfreq(len(t), d=t[1] - t[0])\n",
    "\n",
    "    # Shift the spatial and temporal frequencies to be centered??\n",
    "    spatial_freq_shifted = fftshift(spatial_freq) #don't need it because already have it in : sample_freq_x & sample_freq_shifted_dev\n",
    "    temporal_freq_shifted = fftshift(temporal_freq)\n",
    "\n",
    "    # Plot the magnitude of the 2D Fourier transform\n",
    "    #plt.imshow(np.log(np.abs(c_shifted)), extent=[spatial_freq_shifted[0], spatial_freq_shifted[-1], temporal_freq_shifted[0], temporal_freq_shifted[-1]], aspect ='auto')\n",
    "    #plt.colorbar()\n",
    "    #plt.title('Fourier Transform of contrast')\n",
    "    #plt.ylabel('time freq') # CHECK AGAIN MARIA\n",
    "    #plt.xlabel('spatial freq')\n",
    "    return(c_shifted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### for the f\n",
    "def filtered_contrast_f(sample_freq_x, sample_freq_shifted,c_shifted, f_shifted,h_shifted):\n",
    "    filtered_c_f = np.empty([len(sample_freq_x),len(sample_freq_shifted)],dtype = 'complex_')\n",
    "    for i, fx_value in enumerate(sample_freq_x):# fx\n",
    "        filtered_c_f[i,:] = c_shifted[i,:]*f_shifted[:]*h_shifted[i] #because f_shifted, h_shifted is 1D\n",
    "    # now invert the FT\n",
    "    filtered_c_f = fftshift(filtered_c_f)\n",
    "    filtered_c_f = ifft2(filtered_c_f)\n",
    "    filtered_c_f = filtered_c_f.real\n",
    "    return(filtered_c_f)\n",
    "\n",
    "def filtered_contrast_f_dev(sample_freq_x, sample_freq_shifted_dev, c_shifted,f_dev_shifted,h_shifted):\n",
    "    filtered_c_dev = np.empty([len(sample_freq_x),len(sample_freq_shifted_dev)],dtype = 'complex_')\n",
    "    for i, fx_value in enumerate(sample_freq_x):# fx\n",
    "        filtered_c_dev[i,:] = c_shifted[i,:]*f_dev_shifted[:]*h_shifted[i] \n",
    "    # now invert the FT\n",
    "    filtered_c_dev = fftshift(filtered_c_dev)\n",
    "    filtered_c_dev = ifft2(filtered_c_dev)\n",
    "    filtered_c_dev = filtered_c_dev.real\n",
    "    return(filtered_c_dev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ramp(y):\n",
    "    return np.maximum(0, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c=np.sin()  #g_l=1/10 #g_inh = 0.3*g_l #g_l=1/10 #g_exc = 0.1*g_l\n",
    "def g_vm_cal_calculation(filtered_c,filtered_c_dev):\n",
    "    def ramp(y):\n",
    "        return np.maximum(0, y) # ramp function is contrast dependent\n",
    "    E_exc = 60 # mV\n",
    "    E_inh = -30 # mV\n",
    "    x=np.arange(-180, 180, 0.5) # x=np.arange(-180, 180, 0.5) \n",
    "    V_m = np.empty([len(x),len(t)])\n",
    "    calcium = np.empty([len(x),len(t)])\n",
    "    g1_div_g_l = np.empty([len(x),len(t)])\n",
    "    g2_div_g_l = np.empty([len(x),len(t)])\n",
    "    g3_div_g_l = np.empty([len(x),len(t)])\n",
    "    for i,x_value in enumerate(x[:720]):  #s should be an array calculated over x so the index changes by x (was 710 but works for 720 or len(x))\n",
    "\n",
    "        g1_index = (i - 10) % 720 # Here I make sure that the peripheral neurons get input from spatially seperated points in space \n",
    "        g3_index = (i + 10) % 720 # (index circles around so if i=719 then g3 gets input from i=9 beacuse of this mod )\n",
    "\n",
    "        g1_div_g_l[i, :] = 0.3 * (ramp(-filtered_c[g1_index, :]))\n",
    "        g2_div_g_l[i, :] = 0.1 * (ramp(filtered_c_dev[i, :]))\n",
    "        g3_div_g_l[i, :] = 0.3 * (ramp(filtered_c[g3_index, :]))\n",
    "\n",
    "        # actually i cannot compute the g1,2,3 seperatly because the membrane potential is always g1/g_L\n",
    "        ##g1_div_g_l[i,:] =  0.3*(ramp(-filtered_c[i-10,:]))          #Mi9 g1/g_l   Δx is 5 deg\n",
    "        ##g2_div_g_l[i,:] =  0.1*(ramp(filtered_c_dev[i,:]))          #Tm3 g2/g_l\n",
    "        ##g3_div_g_l[i,:] =  0.3*(ramp(filtered_c[i+10,:]))           #Mi4 g3/g_l\n",
    "        #V_m[i,j] = (g1_div_g_l[i,j]*E_inh+g2_div_g_l[i,j]*E_exc+g3_div_g_l[i,j]*E_inh)/(1+0.3+0.1+0.3) # incorporate my LIF code but first make sure the g are right\n",
    "        \n",
    "        ##### Calculate Membrane Potential\n",
    "        V_m[i,:] = (g1_div_g_l[i,:]*E_inh + g2_div_g_l[i,:]*E_exc + g3_div_g_l[i,:]*E_inh)/(1+g1_div_g_l[i,:]+g2_div_g_l[i,:]+g3_div_g_l[i,:])\n",
    "        calcium[i,:] = (ramp(V_m[i,:]))**2\n",
    "    return(g1_div_g_l,g2_div_g_l,g3_div_g_l,V_m,calcium)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "lif",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}