21102022_FATES_external_dynamic_biases_sensitivity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9eaec9c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import scipy.special as sp\n",
    "import random as rand\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "from numba import jit,int64,float64\n",
    "import time\n",
    "import scipy.optimize as sciopt\n",
    "%matplotlib inline\n",
    "\n",
    "#To increase cell width:\n",
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "\n",
    "#Color-blind friendly colors to plot:\n",
    "CB_colors = ['#377eb8', '#ff7f00', '#4daf4a','#f781bf', '#a65628', '#984ea3','#999999', '#e41a1c', '#dede00']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "627dfa58",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Parameters - self-activation a, cross-repression r, decay k, noise alpha:\n",
    "a = np.array([0,0])\n",
    "r = np.array([1,1])\n",
    "k = np.array([1,1])\n",
    "alpha = np.array([0.1,0.1])\n",
    "#Hill coefficients:\n",
    "n = np.array([[4,4],[4,4]])\n",
    "theta = np.array([[0.5,0.5],[0.5,0.5]])\n",
    "\n",
    "#The total time and dt::\n",
    "T = 40\n",
    "dt = 0.005\n",
    "#Other time-related variables related to these two:\n",
    "Nt = round(T/dt)\n",
    "sqrt_dt = np.sqrt(dt)\n",
    "TimeRange = np.linspace(0,T,num=Nt)\n",
    "# #This is to calculate the dynamic threshold moving average:\n",
    "# timeBase = 2*int(Nt/100)\n",
    "# #Time to plot utility against:\n",
    "# time_trunc = TimeRange[int((timeBase/2)-1):-int(timeBase/2)]\n",
    "\n",
    "#Let's call num_traj the number of trajectories:\n",
    "num_traj = 10000\n",
    "\n",
    "#The threshold above which we say x_i is high:\n",
    "thresh = np.array([0.8,0.8])\n",
    "#The stability threshold; a trajectory has to stay above this:\n",
    "stability_thresh = 0.9\n",
    "\n",
    "#A minimum tolerance:\n",
    "eps = 1e-6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f237f36f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Self-activation:\n",
    "@jit(nopython=True)\n",
    "def activation(x,a,n,theta):\n",
    "    if (x>=0):\n",
    "        return (a*x**n)/(x**n + theta**n)\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "#Cross-inhibition\n",
    "@jit(nopython=True)\n",
    "def repression(x,r,n,theta):\n",
    "    if (x>0):\n",
    "        return (r*theta**n)/(x**n + theta**n)\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "# #Bias (for external signals):\n",
    "# @jit(nopython=True)\n",
    "# def ext_bias(x1,x2,scale1,scale2,t):\n",
    "# #     return (x1-x2,x2-x1)\n",
    "#     temp = 0\n",
    "#     return (temp,0)\n",
    "\n",
    "#Constant external bias:\n",
    "@jit(nopython=True)\n",
    "def ext_bias_const(field):\n",
    "    return (field,0)\n",
    "    \n",
    "# #This solves the system fwd using simple Euler-Maruyama:\n",
    "# @jit(nopython=True)\n",
    "# def Solver(initial,a,r,k,n,theta,alpha):\n",
    "#     final = np.empty((2,num_traj,Nt))\n",
    "#     final[:,:,0] = initial\n",
    "#     #Solving the system forward in time:\n",
    "#     for i in range(num_traj):\n",
    "#         for t in range(1,Nt):\n",
    "#             #Equation for first species:\n",
    "#             temp = np.sqrt(np.maximum(final[0,i,t-1],0.01))\n",
    "#             noise = rand.normalvariate(0,alpha[0]*temp)\n",
    "#             final[0,i,t] = final[0,i,t-1] + dt*(activation(final[0,i,t-1],a[0],n[0,0],theta[0,0]) \\\n",
    "#                                                + repression(final[1,i,t-1],r[0],n[1,0],theta[1,0]) \\\n",
    "#                                                - k[0]*final[0,i,t-1] + ext_bias(final[0,i,t-1],final[1,i,t-1],t)[0]) \\\n",
    "#                                                 + sqrt_dt*noise\n",
    "#             #Equation for second:\n",
    "#             temp = np.sqrt(np.maximum(final[1,i,t-1],0.01))\n",
    "#             noise = rand.normalvariate(0,alpha[1]*temp)\n",
    "#             final[1,i,t] = final[1,i,t-1] + dt*(activation(final[1,i,t-1],a[1],n[1,1],theta[1,1]) \\\n",
    "#                                                + repression(final[0,i,t-1],r[1],n[0,1],theta[0,1]) \\\n",
    "#                                                - k[1]*final[1,i,t-1] + ext_bias(final[0,i,t-1],final[1,i,t-1],t)[1]) \\\n",
    "#                                                 + sqrt_dt*noise\n",
    "#     return final\n",
    "\n",
    "#A solver with some constant extrinsic bias:\n",
    "@jit(nopython=True)\n",
    "def Solver_const_ext_bias(initial,a,r,k,n,theta,alpha,ext_bias_field):\n",
    "    final = np.empty((2,num_traj,Nt))\n",
    "    final[:,:,0] = initial\n",
    "    #Solving the system forward in time:\n",
    "    for i in range(num_traj):\n",
    "        for t in range(1,Nt):\n",
    "            #Equation for first species:\n",
    "            temp = np.sqrt(np.maximum(final[0,i,t-1],0.01))\n",
    "            noise = rand.normalvariate(0,alpha[0]*temp)\n",
    "            final[0,i,t] = final[0,i,t-1] + dt*(activation(final[0,i,t-1],a[0],n[0,0],theta[0,0]) \\\n",
    "                                               + repression(final[1,i,t-1],r[0],n[1,0],theta[1,0]) \\\n",
    "                                               - k[0]*final[0,i,t-1] + ext_bias_const(ext_bias_field)[0]) \\\n",
    "                                                + sqrt_dt*noise\n",
    "            #Equation for second:\n",
    "            temp = np.sqrt(np.maximum(final[1,i,t-1],0.01))\n",
    "            noise = rand.normalvariate(0,alpha[1]*temp)\n",
    "            final[1,i,t] = final[1,i,t-1] + dt*(activation(final[1,i,t-1],a[1],n[1,1],theta[1,1]) \\\n",
    "                                               + repression(final[0,i,t-1],r[1],n[0,1],theta[0,1]) \\\n",
    "                                               - k[1]*final[1,i,t-1] + ext_bias_const(ext_bias_field)[1]) \\\n",
    "                                                + sqrt_dt*noise\n",
    "    return final\n",
    "\n",
    "\n",
    "#Classifier:\n",
    "# @jit(nopython=True)\n",
    "def fate_classifier(traj):\n",
    "    #Stability factor of trajectories:\n",
    "    cross_flags = np.zeros((2,num_traj))\n",
    "    cross_times = np.ones((2,num_traj))*Nt\n",
    "    for axis_idx in range(2):\n",
    "        #Axis crossings:\n",
    "        for traj_idx in range(num_traj):\n",
    "            if (np.size(np.where(traj[axis_idx,traj_idx]>thresh[axis_idx])[0]) != 0):\n",
    "                cross_flags[axis_idx,traj_idx] = 1\n",
    "                cross_times[axis_idx,traj_idx] = np.where(traj[axis_idx,traj_idx]>thresh[axis_idx])[0][0]\n",
    "                    \n",
    "    #Stability factor: after the threshold is crossed, how much time the traj spends above it:\n",
    "    stability_factrs = np.zeros((2,num_traj))\n",
    "    for axis_idx in range(2):\n",
    "        for traj_idx in range(num_traj):\n",
    "#             if (cross_flags[axis_idx,traj_idx]==0):\n",
    "#                 stability_factrs[axis_idx,traj_idx] = -1\n",
    "            if (cross_flags[axis_idx,traj_idx]==1):\n",
    "                stability_factrs[axis_idx,traj_idx] = np.sum(traj[axis_idx,traj_idx,int(cross_times[axis_idx,traj_idx]):]>thresh[axis_idx])\\\n",
    "                /len(traj[axis_idx,traj_idx,int(cross_times[axis_idx,traj_idx]):])\n",
    "    \n",
    "    #Stability threshold - trajectories that spend more than this above the concentration threshold are considered committed:\n",
    "#     stability_thresh = 0.8\n",
    "    #Classifying fates:\n",
    "    fates = np.zeros((2,num_traj))\n",
    "    for axis_idx in range(2):\n",
    "        for traj_idx in range(num_traj):\n",
    "            if (cross_times[axis_idx,traj_idx]<=int(Nt/2) and stability_factrs[axis_idx,traj_idx]>=stability_thresh):\n",
    "                fates[axis_idx,traj_idx]=1\n",
    "                \n",
    "    return stability_factrs,fates\n",
    "\n",
    "def fate_fractions(fates):\n",
    "    #Initializing:\n",
    "    fate_frax = np.zeros(4)\n",
    "    fate_frax[0] = np.sum((fates[0]==0) & (fates[1]==0))/num_traj\n",
    "    fate_frax[1] = np.sum((fates[0]==1) & (fates[1]==0))/num_traj\n",
    "    fate_frax[2] = np.sum((fates[0]==0) & (fates[1]==1))/num_traj\n",
    "    fate_frax[3] = np.sum((fates[0]==1) & (fates[1]==1))/num_traj\n",
    "    \n",
    "    return fate_frax\n",
    "\n",
    "def traj_moments(traj,fates):\n",
    "    #Finding the fate fractions first:\n",
    "    fate_frax = fate_fractions(fates)\n",
    "    \n",
    "    #Flags, if fate_frac is zero for a fate then no avg or std is calculated:\n",
    "    fate_flags = np.array([False for i in range(4)])\n",
    "    \n",
    "    #Flag = True if that fate exists in the population:\n",
    "    for fate_idx in range(4):\n",
    "        if (fate_frax[fate_idx] != 0):\n",
    "            fate_flags[fate_idx] = True\n",
    "            \n",
    "    #Average trajectories and standard deviation around them:\n",
    "    avg_traj = np.zeros((4,2,Nt))\n",
    "    std_traj = np.zeros((4,2,Nt))\n",
    "    \n",
    "    if (fate_flags[0]==True):\n",
    "        avg_traj[0,:,:] = np.mean(traj[:,(fates[0]==0) & (fates[1]==0),:],axis=1)\n",
    "        std_traj[0,:,:] = np.std(traj[:,(fates[0]==0) & (fates[1]==0),:],axis=1)\n",
    "\n",
    "    if (fate_flags[1]==True):\n",
    "        avg_traj[1,:,:] = np.mean(traj[:,(fates[0]==1) & (fates[1]==0),:],axis=1)\n",
    "        std_traj[1,:,:] = np.std(traj[:,(fates[0]==1) & (fates[1]==0),:],axis=1)\n",
    "\n",
    "    if (fate_flags[2]==True):\n",
    "        avg_traj[2,:,:] = np.mean(traj[:,(fates[0]==0) & (fates[1]==1),:],axis=1)\n",
    "        std_traj[2,:,:] = np.std(traj[:,(fates[0]==0) & (fates[1]==1),:],axis=1)\n",
    "\n",
    "    if (fate_flags[3]==True):\n",
    "        avg_traj[3,:,:] = np.mean(traj[:,(fates[0]==1) & (fates[1]==1),:],axis=1)\n",
    "        std_traj[3,:,:] = np.std(traj[:,(fates[0]==1) & (fates[1]==1),:],axis=1)\n",
    "        \n",
    "    return (avg_traj,std_traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1950d6f2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEjCAYAAADkAazgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhOklEQVR4nO3de5wcVZ338c/XAMqicjFcE+IEDI9Edg1hiIgidzdcJIruIygQcCWiRlFZNQgq4D4ruiiCcjEsSEAkIhCMPpEAAbwhSMIdAhJDgJBAgtxvhsBv/zinTafTPdPTUzU9M3zfr1e9evqcU1W/rump39SpU1WKCMzMzMr0unYHYGZmg5+TjZmZlc7JxszMSudkY2ZmpXOyMTOz0jnZmJlZ6ZxszHpB0m6SQtIJTbY/PLc/vAfrOD/P09FimK85Pf29WPmcbOwf8h9nd9NuvVi+d5pWCEkd+bt0frtjseas1e4ArF86sYu6RX0VhP3DscDJwCPtDmQA+TOwLfB4uwOxxMnG1hARJ7Q7BlslIpYCS9sdx0ASES8A97Y7DlvF3WjWEkkjJT0l6QlJb62pW0/SfEmvSNo1lwUwMTd5oKpbblHNvBtJ+nae/0VJT0uaI+n9dWL4x/kPSeMlXZ/bR5363XP9s5KekfT/JW1bZ5nbSDpZ0lxJyyX9XdKDkqZKGl7Q5qusaz9JN0h6XtKTki6VNKpOu7rdj/lzXSZpYd5Wz0j6o6RDGqxvq/w5FuT2T0i6U9LZkt7STazD8u/zli7aXJnj3K6q7ID8+1uat+USSb+V9JluN1Dj9ZwAPJDfTqzp5j08t6l7ziZ/B0LS2pK+Iemvkl6SdK+kI6vaHZW3zYuSFks6UVLd/aWkd+Xf3aOSVkh6WNKPJW3R6mccjHxkYy2JiAckfRL4BXCxpPdFxMpcfSbwduCEiPhtLjsR+CDwTuA04KlcXnklJ63rgQ7g98CVwHrA/sCVkj4VEefUCecjwHjgN8DZef5q+wMTqupHA/sCO0oaHRHVXS0HAkcB1wE3ACuAdwCfBD4gqTMiiujOOhDYB5hB+sxjgA8Du0vaOSLua2IZZwH3AL8jHfm8JX+uCyX9n4j4eqWhpM2Bm4E3A7OAy4A3ACOBQ4EfAX9rtKKIeETSNcD7Jf1zRNxZXZ+XvxcwLyLuymWTgB8DjwK/InVpbQL8C3AE6XvSiuuBDYCjgduBK6rqbmtyGdOBd5G2xcuk79BUSS/n+CYCvwbmAAcA3wBeAL5TvRBJRwDnAH8HZgIPA6NY9X3ZKSIe6tnHG6QiwpMnIgIg8nRCg2lKnXnOzPN8O78/LL+/DnhdTdvzc11Hg/VfD7wKHFRTvgFpJ/IisGlV+eF5ea8C4+ssr1K/Etizpu7bue4rNeXDgNfXWdb7gVeAs2rKd6tssya3cSWmAPavqTs6l89pZrsBW9dZ/jqkHeTLwLCq8s/lZRxdZ571gHWbiP3gvIxT6tR9Odd9rqpsHmknvEmd9kN7+V3tyOs7v0F93d9L/o4FKfFuUFW+FekfiydJR03V224DUqJcDqxVVb5NnmdBdftct0f+vszozeccTFPbA/DUf6aqnWCj6ak687whJ4JXgcnAc8AyYIs6bevuNHPdO3PdLxrENiHXf6aqrLLjntFgnkr9T+vUjcx1l/Zg+9wBLKwpq7tT62IZlZjm1KkbkndcAby1me3WYB0H5vaHVZVVks2kXnw/1iUdiS4FhtTU3ZV3vEOryuYBzwMblvBd7W2y2bPOPNfmuk/UqftJnd/LqblsvwYxzCD9o/Omoj//QJzcjWZriAj1oO1Lkj4KzAV+SPrj+0hELOnhat+dX9ev7WfPNs6va5xnIY086srcOmUP59cNqwslCfg4KSm8M9cPqWqyopt1Neu3tQUR8YqkPwBbA9sDD3a1AEkjgK8CewIjSMmg2rCqn2cC/wWcIelfgdnAH4F7Iu8ZuxMRL0q6BDgS+FdSFxSSdiB1Nc6I1bskLwK+B9wt6ef5M/8xIpY3s76S1ftOVL6z8+rUVbpOh7Pq91L5zu4qacc682xC+u5s02CZrylONlaEv5D+69+ZdA7hqhaWUTlBvXeeGnljnbJHu1n2U7UFEbEy5ZXVEgnA94EvkP57n03aybyY6w4H3koxHmtQXvks63c1s6StSEl2Q9L5rauAp0ldNx2kcw6vr7SPiAcljSN1h44nHf0APCzplIg4vcm4zyclm4nkZMOqgR/TqhtGxPclPQ58Bvg8abuGpN8CX46Iejv8PhERT9cprpxz7Kpu7aqyynf2y92srt539jXHycaKMIWUaB4n/Yd7LPD/eriMyh/40T3Y8VU09Z95dyRtQtop3gXsHBHP1tQfXMR6sk0blG+WX+vt8Kp9ibSzOyIizq+uyHFOrJ0hIuYDH5W0FumobS9S99ppkp6PiHO7CzoibpB0PzBB0gakbrKDSb/7WXXaXwBckNvuDHwI+AQwW9K2EbGsu3X2Y5Xf0foR8UxbIxkAPPTZekXSzsBJwH3Advn1REnvrdP8lfxaezQBcGN+3aXwIJu3Felv4qo6iWZ4ri/KrrUFkoYAle12azfzvy2/XtbMsqtFxMqImBcR3yElCkgjBZs1jXTU9FFgP2Ao8LOIeLmLdT4VEbMi4kjS0dFG9O533dV3qa/0h+/sgOFkYy2TtCFwMekP/6CIeIy0A1pJGg5de+1GZWjtiNpl5S6V3wMHSvpEg/X9cz76KMui/PrevOOvrPeNpOGtRfYE7CFp/5qyyaTzNddFRJfna1gV627Vhfl8zCdrG0saJ6ne0VSl7IXuAq5yAWlAyGF5gpRAatc5Ph9F1ar8Dl+oajtU0tslDW0yhidJR7RrfJf60I9Io/5OlbRNbaWkdSQ5EWXuRrM1NDhBX3FFRNyWfz6P9Mf++UpZRNwu6RjSH+JPSNcoVMwh9W+fI+lS0si1pyLiR7n+Y6QRQedK+jxwE+l8y3DStQ/bkU7KltL1EhGPSpoOHATcJukq0rmTvYGXSKPuxhS0ul8BMyTNII1AeyfpGpknSOc4unMm6VqVX0i6jHRuaTvS+ZhLSEm/2seAz+bzJQtIO+utgQ+Qhif/oNnAI+JhSdeRBiasBO6MiHpHYtOBl/Kgh0WASEcBO5JOmF9T1XYy8E3S9VgnNBHDc5JuAnaRdBHpvOErwMyIuKPZz9IbEXFv/sfoPNIgiCtzHGuT/i52IQ2XfntfxNPfOdlYPd/som4RaUf8OVLXy8yI+GF1g4g4Q9KewIckfTEiTs3ls3MiOhL4IumakAdJiYmIWJxHNn2OdIHjx0ndJI+SBh78EFjtYsIS/DuwkLSz/ixpZzGTdFFfvS6rVl0OTAWOI3VFvZzLjo2Iv3Q3c0TcIWl34D9JSWot0gWOB5ISdG2yuZjU9bUzMJY0cu0RUkL4XuQLMXvgfFKyWYuagQFVppBGrY3NMb5E+n1/lXS9UsNutyYdShp+PJ7UHShgMWmwSp+IiJ9Kuh04BtiddD3W86SRbZcCP++rWPo7NTnq0czMrGU+Z2NmZqVzsjEzs9I52ZiZWemcbMzMrHQejdbA0KFDo6Ojo91hmJkNKPPmzXs8IjauLXeyaaCjo4O5c9t26yYzswFJUt0Lkt2NZmZmpXOyMTOz0jnZmJlZ6ZxszMysdE42ZmZWOicbMzMrnZONmZmVzsnGzMxK52RjZmalc7IxM7PSOdmYmVnpnGzMzKx0TjZmZlY6JxszMyudk42ZmZXOycbMzErnZGNmZqVzsjEzs9I52ZiZWemcbMzMrHRONmZmVjonGzMzK52TjZmZlc7JxszMSudkY2ZmpRswyUbSeEn3SVogaUqdekk6PdffIWlsTf0QSbdK+nXfRW1mZjBAko2kIcAZwD7AaOBgSaNrmu0DjMrTJOCsmvqjgfklh2pmZnUMiGQDjAMWRMTCiFgBTAcm1LSZAFwQyY3ABpI2B5A0HNgP+J++DNrMzJKBkmyGAQ9XvV+cy5pt8wPgK8CrXa1E0iRJcyXNXb58ea8CNjOzVQZKslGdsmimjaT9gWURMa+7lUTE1IjojIjOjTfeuJU4zcysjoGSbBYDW1a9Hw4sabLNe4ADJC0idb/tIemn5YVqZma1BkqyuRkYJWmkpHWAg4CZNW1mAoflUWk7AU9HxNKIODYihkdER57v2og4pE+jNzN7jVur3QE0IyJWSpoMzAaGAOdFxN2Sjsr1ZwOzgH2BBcALwBHtitfMzFaniNpTHwbQ2dkZc+fObXcYZmYDiqR5EdFZWz5QutHMzGwAc7IxM7PSOdmYmVnpnGzMzKx0TjZmZlY6JxszMyudk42ZmZXOycbMzErnZGNmZqVzsjEzs9I52ZiZWemcbMzMrHRONmZmVjonGzMzK52TjZmZlc7JxszMSudkY2ZmpXOyMTOz0jnZmJlZ6ZxszMysdE42ZmZWOicbMzMrnZONmZmVrrBkI2nbopZlZmaDS5FHNtdJ2rG2UNKmBa7DzMwGoCKTzRHAryXtBSBpPUknAvcXuA4zMxuA1ipqQRHxG0kfBC6V9HPgY8B8YHxR6zAzs4Gp6AECw4CXgKOB6RGxe0TcUPA6zMxsgCnsyEbSrcBGwDeA24BfSbonIqYWtQ4zMxuYCks2wMXAaRHxdwBJuwKzJW0aEd8qcD1mZjbAFNaNFhHfrSSa/P4BYBfgQ0UsX9J4SfdJWiBpSp16STo9198haWwu31LSdZLmS7pb0tFFxGNmZs0r9aLOiHgM2LW3y5E0BDgD2AcYDRwsaXRNs32AUXmaBJyVy1cCx0TEtsBOwGfrzGtmZiVqKdlIGiJpa0lbdNc2Ip5tZR01xgELImJhRKwApgMTatpMAC6I5EZgA0mbR8TSiLilKpb5pIEMZmbWR3p8zkbSe4BLgU3y+2eBW4FbqqZ7IyIKjHMY8HDV+8XAu5poMwxYWhV7B7A9cFOBsZmZWTdaGSBwGvBG4L+AvwPbATsA7wMEBPA88OaCYiQvt1ZtMuuyjaQ3ApcBX4iIZ+quRJpE6oJjxIgRrUVqZmZraCXZbAt8JyJOqi6UtD4wlpR4xhYQW7XFwJZV74cDS5ptI2ltUqK5KCIub7SSPEx7KkBnZ2eRR2ZmZq9prSSbR4HltYUR8TRwXZ6KdjMwStJI4BHgINIdCqrNBCZLmk7qYns6IpZKEnAuMD8ivl9CbGZm1o1WBgj8HNit4Di6FBErgcnAbNIJ/ksi4m5JR0k6KjebBSwEFgDnAJ/J5e8BDgX2kHRbnvbty/jNzF7r1N15fEl3kO4IcEeeHgB+AXw3In5WdoDt0tnZGXPnzm13GGZmA4qkeRHRWVveTDfay8C/AYew6oT734ELJX2YdKRzU0Q8WFSwZmY2uHSbbCJih3xR5WhgDGno8Jg8fShPIelJVg2BnhcRl5QTspmZDTRNDRCIiFeAO/N0YaU8X7cyhlUJaHtgT9IRkJONmZkBvbwRZ0QsAhYBV1TKJG1ESjpmZmZAsXd9BiAingDmFL1cMzMbuEq9EaeZmRk42ZiZWR9wsjEzs9I52ZiZWemcbMzMrHRONmZmVrpuhz5LuraF5UZE7NnCfGZmNgg1c53NVqz5oLL1gKH556dIDy5bP79/HHiuiODMzGxw6LYbLSI6ImJkZSLdjuZF0hM7t4iIjSJiQ2AL4HTghdzGzMwMaO2czanADRHxxYh4tFIYEY9GxBeAG3MbMzMzoLVksxvw2y7qrwd2byUYMzMbnFpJNgFs20X9O1jzHI+Zmb2GtZJsrgI+LekwSaoUKpkIfCq3MTMzA1q76/OXgB2BnwAnS7qfdCSzDbAp8HBuY2ZmBrRwZBMRi0kPSvsO8CQwDnhX/vk7wJjcxszMDGjxeTYR8TTwtTyZmZl1yberMTOz0jVzu5rDWllwRFzQynxmZjb4NNONdj5pAIC6aVctACcbMzMDmks2vkDTzMx6pdtkExFd3S3AzMysW70aICDp9ZKGSVqnqIDMzGzwaSnZSBqbn3PzLPAQ8N5cvomkOZL2KjBGMzMb4HqcbCSNAX4PbE3NIICIWAasC0wsIjgzMxscWjmyOQlYQrrh5hTWHKU2h3RXATMzM6C1ZLMLcE5EPEf9uzs/RHqQmpmZGdBasnkD8HQX9W9uMZYuSRov6T5JCyRNqVMvSafn+jskjW12XjMzK1cryeavwA5d1O8B3NNaOPVJGgKcAewDjAYOljS6ptk+wKg8TQLO6sG8ZmZWolZuxPkz4OuSLgFuzWUBIOkYYDxwdDHh/cM4YEFELMzrmQ5MYPWkNgG4ICICuFHSBpI2BzqamLcw55xzDrNnzy5j0WZmfeK4445j++23L3SZrSSbU4C9gdnAvaREc6qkjYHNgKuBMwuLMBlGek5OxWLSYw26azOsyXkBkDSJdFTEiBEjWgr0scce4957721pXjOz/uD5558vfJk9TjYRsULS3sDngI8DL5EenHY/8H3gtIh4tdAo69+XrXZwQqM2zcybCiOmAlMBOjs7W3q09fHHH8/xxx/fyqxmZoNWq8+zWQmcmqe+sBjYsur9cNLw62barNPEvGZmVqKB8jybm4FRkkbmW+McBMysaTMTOCyPStsJeDoiljY5r5mZlagnz7O5MCKi2efbFPk8m4hYKWky6TzREOC8iLhb0lG5/mxgFrAvsAB4ATiiq3mLis3MzLqnNHiriwbSq6RzHOvm8zWV91093yYiYkhxYfa9zs7OmDt3brvDMDMbUCTNi4jO2vKmn2cTESuq35uZmTWrmWQzEfhx1fsA5kfE8nJCMjOzwaaZAQKHk+7wXHEd6TobMzOzpjSTbB4HNq1639W5GjMzszU00412A3C8pBHAk7nsQElv62KeiIhv9To6MzMbFJpJNl8ApgGfJx3VBHBgnhoJwMnGzMyAJpJNRCwCds0XRG4GLCIloF+WGZiZmQ0eTd+uJg99fkjSNOCmiHiwvLDMzGwwaeVGnEeUEYiZmQ1eA+XeaGZmNoA52ZiZWemcbMzMrHRONmZmVjonGzMzK52TjZmZlc7JxszMSldospG0paRrJc0pcrlmZjaw9fiizm78E7Ab6d5oZmZmQPHJ5q/AyIKXaWZmA1yhySYiVgK+Z5qZma2mx+dsJJ0n6V1d1I+TdF7vwjIzs8GklQECh7P6Y6JrjQQmthSNmZkNSmUMfV4PeLmE5ZqZ2QDV1Dmb/Ejojqqit0t6X52mGwGfBhb0PjQzMxssmh0gcATwTdKQ5gCOy1MtAa/m9mZmZkDzyeYK0uOgBZwHTAX+VNMmgOeAmyPi4YLiMzOzQaCpZBMRtwO3A0h6K3B5RNxZZmBmZjZ4tDJAQMBdDSuljSRd0XJEZmY26LSSbL4BXCdpWG2FpF1JR0D79DYwMzMbPFpJNkcBOwK3S5oAIOl1kr4FXAOsBHYtLkQzMxvoepxsImIqMA5YClwu6Wzgd6TRaZcDYyLixqICzN1yV0u6P79u2KDdeEn3SVogaUpV+X9LulfSHZJmSNqgqNjMzKw5LV3UGRF3k45u/gQcCbwb+FpEfDQini4wPoApwJyIGAXMye9XI2kIcAap+240cLCk0bn6amC7iPgX4C/AsQXHZ2Zm3Wgp2UhaB/gusDOwkNR1NjmfsynaBGBa/nka8ME6bcYBCyJiYUSsAKbn+YiIq/INQgFuBIaXEKOZmXWhlRtxbkPaaU8GzgK2A3YBVgDXSDpJUpG3wdk0IpYC5NdN6rQZBlRf27M4l9X6BPCbAmMzM7MmtPKIgXmke599OCJm5LI/S3on6WLP40kPUKt3O5u6JF0DbFanqt5dCuouok7Zag9wk3Qc6Qjsoi7imARMAhgxYkSTqzYzs+60kmxuBz4WEQ9VF0bEs6RzJdcAp/VkgRGxV6M6SY9J2jwilkraHFhWp9liYMuq98OBJVXLmAjsD+wZEQ2fIpoHP0wF6Ozs9NNGzcwK0kp31661iaZaRJwL7NB6SGuYyapHFkwEflmnzc3AKEkj8/mkg/J8SBoPfBU4ICJeKDAuMzNrUitDn1+p/Czp9ZKG5R18dZv7igguOxnYW9L9wN75PZK2kDQrr28l6RzSbGA+cEkeMQfwI+BNwNWSbstDtc3MrA+19FhoSWOBU4D3AkNISeBaSZsAFwPfjohriggwIv4G7FmnfAmwb9X7WcCsOu3eVkQcZmbWulZGo40Bfk96WucF1XURsQxYFz+p08zMqrRyzuYk0sn3d5AusKwdCTaHdN2LmZkZ0Fqy2QU4JyKeo2Z4cfYQsEWvojIzs0GllWTzBqCrW9K8ucVYzMxskGol2fyVroc27wHc01o4ZmY2GLWSbH4GHCqp+kLMAJB0DDAeuLCA2MzMbJBoZejzKaShzrOBe0mJ5lRJG5NuOXM1cGZhEZqZ2YDXykWdK0jJ5j+AF4GXgG2Ax4GvAPtHxKtFBmlmZgNbSxd15iv2T82TmZlZl4p8FICZmVldTjZmZlY6JxszMyudk42ZmZXOycbMzErnZGNmZqVzsjEzs9I52ZiZWemcbMzMrHRONmZmVjonGzMzK52TjZmZlc7JxszMSudkY2ZmpXOyMTOz0jnZmJlZ6ZxszMysdE42ZmZWOicbMzMrnZONmZmVzsnGzMxK52RjZmal6/fJRtJGkq6WdH9+3bBBu/GS7pO0QNKUOvX/ISkkDS0/ajMzq9bvkw0wBZgTEaOAOfn9aiQNAc4A9gFGAwdLGl1VvyWwN/BQn0RsZmarGQjJZgIwLf88DfhgnTbjgAURsTAiVgDT83wVpwJfAaLEOM3MrIGBkGw2jYilAPl1kzpthgEPV71fnMuQdADwSETc3t2KJE2SNFfS3OXLl/c+cjMzA2CtdgcAIOkaYLM6Vcc1u4g6ZSHpn/Iy3t/MQiJiKjAVoLOz00dBZmYF6RfJJiL2alQn6TFJm0fEUkmbA8vqNFsMbFn1fjiwBNgaGAncLqlSfoukcRHxaGEfwMzMujQQutFmAhPzzxOBX9ZpczMwStJISesABwEzI+LOiNgkIjoiooOUlMY60ZiZ9a2BkGxOBvaWdD9pRNnJAJK2kDQLICJWApOB2cB84JKIuLtN8ZqZWY1+0Y3WlYj4G7BnnfIlwL5V72cBs7pZVkfR8ZmZWfcGwpGNmZkNcE42ZmZWOicbMzMrnZONmZmVzsnGzMxK52RjZmalc7IxM7PSOdmYmVnpnGzMzKx0TjZmZlY6JxszMyudk42ZmZXOycbMzErnZGNmZqVzsjEzs9I52ZiZWemcbMzMrHRONmZmVjonGzMzK52TjZmZlc7JxszMSudkY2ZmpXOyMTOz0jnZmJlZ6RQR7Y6hX5K0HHiwxdmHAo8XGE5RHFfPOK6ecVw901/jgt7F9taI2Li20MmmBJLmRkRnu+Oo5bh6xnH1jOPqmf4aF5QTm7vRzMysdE42ZmZWOiebckxtdwANOK6ecVw947h6pr/GBSXE5nM2ZmZWOh/ZmJlZ6ZxszMysdE42BZM0XtJ9khZImtLueCokLZJ0p6TbJM1tYxznSVom6a6qso0kXS3p/vy6YT+J6wRJj+RtdpukfdsQ15aSrpM0X9Ldko7O5W3dZl3E1dZtJukNkv4s6fYc14m5vN3bq1Fcbf+O5TiGSLpV0q/z+8K3l8/ZFEjSEOAvwN7AYuBm4OCIuKetgZGSDdAZEW29iEzS+4DngAsiYrtc9l3giYg4OSfoDSPiq/0grhOA5yLilL6MpSauzYHNI+IWSW8C5gEfBA6njdusi7j+L23cZpIErBcRz0laG/gDcDRwIO3dXo3iGk+bv2M5vi8BncCbI2L/Mv4mfWRTrHHAgohYGBErgOnAhDbH1K9ExO+AJ2qKJwDT8s/TSDutPtUgrraLiKURcUv++VlgPjCMNm+zLuJqq0iey2/XzlPQ/u3VKK62kzQc2A/4n6riwreXk02xhgEPV71fTD/4A8wCuErSPEmT2h1MjU0jYimknRiwSZvjqTZZ0h25m63Pu/eqSeoAtgduoh9ts5q4oM3bLHcJ3QYsA66OiH6xvRrEBe3/jv0A+ArwalVZ4dvLyaZYqlPWL/57Ad4TEWOBfYDP5m4j69pZwNbAGGAp8L12BSLpjcBlwBci4pl2xVGrTlxt32YR8UpEjAGGA+MkbdfXMdTTIK62bi9J+wPLImJe2etysinWYmDLqvfDgSVtimU1EbEkvy4DZpC6/PqLx/I5gMq5gGVtjgeAiHgs7yBeBc6hTdss9/FfBlwUEZfn4rZvs3px9ZdtlmN5CriedF6k7durXlz9YHu9Bzggn9OdDuwh6aeUsL2cbIp1MzBK0khJ6wAHATPbHBOS1ssncZG0HvB+4K6u5+pTM4GJ+eeJwC/bGMs/VP7Ysg/Rhm2WTyyfC8yPiO9XVbV1mzWKq93bTNLGkjbIP68L7AXcS/u3V9242r29IuLYiBgeER2k/dW1EXEIJWyvtXq7AFslIlZKmgzMBoYA50XE3W0OC2BTYEbaP7AW8LOIuLIdgUi6GNgNGCppMfBN4GTgEkn/DjwE/Fs/iWs3SWNIXaGLgE/1dVyk/zwPBe7M/f0AX6P926xRXAe3eZttDkzLI0NfB1wSEb+W9Cfau70axXVhP/iO1VP498tDn83MrHTuRjMzs9I52ZiZWemcbMzMrHRONmZmVjonGzMzK52TjVkfkbSbpJB0eLtjMetrTjZmvZCTR1MT0NHueJsh6QhJn293HDa4+Dobs16QdEhN0S7AJNIz3H9fUzcDeAV4OSJe6YPwWiJpOXBjRHyg3bHY4OE7CJj1QkT8tPq9pLVIyeZPtXUDgaS3AUOBG9sdiw0u7kYz6yP1ztlIOjyX7SnpG5IelPSipJsk7ZTb7CrpD5Kel7RU0tcbLP/1kr6m9CTIlyQ9JelXkrZvMr4rgPvz2/+s6gL8Vu8+uZmPbMz6i5NJ99M7DVgHOAaYLWki6YaXU4GLSE/CPEnSA9VHTvkOzFcCOwMXAj8C1geOBP4o6X0R0d3jwKfmGPYHPk16cinAnwr5hPaa5mRj1j8MAXbKT3hF0j2kO+1eCrw7Im7O5ecCDwKfBaq76SaTbiQ6PiJmVwolnUm6k/Apub6hiJgl6VPA8og4u5iPZZa4G82sfzirkmiyyuCCGyuJBiC3+TMwqmb+Q0i30p8naWhlIh0lXQ28N9/avjtjgVtb/RBmjfjIxqx/WFj9JiKezI+EeKBO2yeBt9SUbQusCyzvYh1DWf2x5avJyWk4qx8xmRXCycasf2g0FLrZIdIC7gS+1EWbrhIRwA759ZYm12nWNCcbs8HhfmBj0pMWX21xGZVRa042VjifszEbHC4ANqPBkY2kTZtYxlb59aGigjKr8JGN2eBwGrA38N+S9gCuBZ4BRgB7Ai8Bu3ezjMp5o9PzY5RfIT1C3LcZsV5zsjEbBCLiZUn7AZ8BDgVOzFVLSKPXpjWxmNOB0cBHgKOAhyLiohLCtdcg3xvNzMxK53M2ZmZWOicbMzMrnZONmZmVzsnGzMxK52RjZmalc7IxM7PSOdmYmVnpnGzMzKx0TjZmZla6/wVxF6f3/XvioAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#The initial conditions:\n",
    "x0 = np.zeros((2,num_traj))\n",
    "#Starting at the saddle point (0.5,0.5):\n",
    "# x0 = 0.5*np.ones((2,num_traj))\n",
    "# #Uniform around mean:\n",
    "# init_bias = 0\n",
    "# x0[0,:] = init_bias*np.ones(num_traj)\n",
    "\n",
    "#For a total (constant) bump of 0.1, we need to apply 0.1/T at each timestep:\n",
    "ext_bias_field = 0/T\n",
    "ext_bias_dyn = ext_bias_field*np.ones_like(TimeRange)\n",
    "\n",
    "plt.plot(TimeRange,ext_bias_dyn,color='k')\n",
    "plt.xlabel(r\"Time $t$\",fontsize=18)\n",
    "plt.ylabel(r\"ext. field $h_x$\",fontsize=18)\n",
    "plt.title(\"External bias vs. time\",fontsize=20)\n",
    "\n",
    "print(np.trapz(ext_bias_dyn,x=TimeRange))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bf92316a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time taken: 14.094251871109009 s.\n"
     ]
    }
   ],
   "source": [
    "#Generating some trajectories:\n",
    "start_time = time.time()\n",
    "traj_simple = np.empty((2,num_traj,Nt))\n",
    "traj_simple = Solver_const_ext_bias(x0,a,r,k,n,theta,alpha,ext_bias_field)\n",
    "end_time = time.time()\n",
    "loop_time = end_time - start_time\n",
    "print(f\"Time taken: {loop_time} s.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "41b54519",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.99985846 0.         0.         ... 0.99082691 0.         0.        ]\n",
      " [0.         0.99804305 0.99466046 ... 0.         0.99822671 0.99904567]]\n",
      "(4, 2, 8000)\n",
      "[0.     0.4999 0.5001 0.    ]\n"
     ]
    }
   ],
   "source": [
    "#Now testing the fates and stability factors:\n",
    "stability_factrs,fates = fate_classifier(traj_simple)\n",
    "fate_frax = fate_fractions(fates)\n",
    "avg_traj,std_traj = traj_moments(traj_simple,fates)\n",
    "#Fate fractions:\n",
    "print(stability_factrs)\n",
    "print(np.shape(avg_traj))\n",
    "print(fate_frax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0b663d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Defining our grid points and making the grid:\n",
    "x = np.arange(0, 2*(a[0]+r[0]), 0.1)\n",
    "y = np.arange(0, 2*(a[1]+r[1]), 0.1)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "\n",
    "#The velocity field:\n",
    "x_dot = (a[0]*X**n[0,0])/(X**n[0,0] + theta[0,0]**n[0,0]) + (r[0]*theta[1,0]**n[1,0])/(Y**n[1,0] + theta[1,0]**n[1,0]) - k[0]*X + ext_bias_field\n",
    "y_dot = (a[1]*Y**n[1,1])/(Y**n[1,1] + theta[1,1]**n[1,1]) + (r[1]*theta[0,1]**n[0,1])/(X**n[0,1] + theta[0,1]**n[0,1]) - k[1]*Y\n",
    "\n",
    "# #Only the additional x-bias:\n",
    "# x_dot = 0.01\n",
    "# y_dot = 0\n",
    " \n",
    "# creating plot\n",
    "fig, ax = plt.subplots(figsize =(6.4, 6.4))\n",
    "ax.quiver(X, Y, x_dot, y_dot)\n",
    "#Diagonal:\n",
    "ax.plot([0, 1],[0, 1],transform=ax.transAxes,color='k',ls=':')\n",
    "#Thresholds:\n",
    "ax.axvline(thresh[0],color='k',ls=':')\n",
    "ax.axhline(thresh[1],color='k',ls=':')\n",
    "\n",
    "#Total state space visited, all:\n",
    "ax.scatter(traj_simple[0,(fates[0]==0)&(fates[1]==0),:],traj_simple[1,(fates[0]==0)&(fates[1]==0),:],color=CB_colors[0],marker=\"o\",alpha=0.01)\n",
    "ax.scatter(traj_simple[0,(fates[0]==1)&(fates[1]==0),:],traj_simple[1,(fates[0]==1)&(fates[1]==0),:],color=CB_colors[1],marker=\"o\",alpha=0.01)\n",
    "ax.scatter(traj_simple[0,(fates[0]==0)&(fates[1]==1),:],traj_simple[1,(fates[0]==0)&(fates[1]==1),:],color=CB_colors[2],marker=\"o\",alpha=0.01)\n",
    "ax.scatter(traj_simple[0,(fates[0]==1)&(fates[1]==1),:],traj_simple[1,(fates[0]==1)&(fates[1]==1),:],color=CB_colors[3],marker=\"o\",alpha=0.01)\n",
    "\n",
    "# # Final points:\n",
    "# ax.scatter(traj_simple[0,(fates[0]==0)&(fates[1]==0),-1],traj_simple[1,(fates[0]==0)&(fates[1]==0),-1],color=CB_colors[0],marker=\"o\",alpha=1)\n",
    "# ax.scatter(traj_simple[0,(fates[0]==1)&(fates[1]==0),-1],traj_simple[1,(fates[0]==1)&(fates[1]==0),-1],color=CB_colors[1],marker=\"o\",alpha=1)\n",
    "# ax.scatter(traj_simple[0,(fates[0]==0)&(fates[1]==1),-1],traj_simple[1,(fates[0]==0)&(fates[1]==1),-1],color=CB_colors[2],marker=\"o\",alpha=1)\n",
    "# ax.scatter(traj_simple[0,(fates[0]==1)&(fates[1]==1),-1],traj_simple[1,(fates[0]==1)&(fates[1]==1),-1],color=CB_colors[3],marker=\"o\",alpha=1) \n",
    "\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlabel(r\"Marker (coordinate) $x$\",fontsize=18)\n",
    "ax.set_ylabel(r\"Marker (coordinate) $y$\",fontsize=18)\n",
    "ax.set_title(\"Flow-field and explored region\",fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df7fa95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's see how the trajectories look:\n",
    "fig,ax = plt.subplots(1,2,figsize=[9.6,4.8],sharey=True)\n",
    "\n",
    "ax[0].plot(TimeRange,avg_traj[1,0],c=CB_colors[1],label=\"x\")\n",
    "ax[0].fill_between(TimeRange,avg_traj[1,0]-std_traj[1,0],avg_traj[1,0]+std_traj[1,0],color=CB_colors[1],alpha=0.5)\n",
    "ax[0].plot(TimeRange,avg_traj[1,1],c=CB_colors[2],label=\"y\")\n",
    "ax[0].fill_between(TimeRange,avg_traj[1,1]-std_traj[1,1],avg_traj[1,1]+std_traj[1,1],color=CB_colors[2],alpha=0.5)\n",
    "ax[0].axhline(thresh[0],color='k',ls='--')\n",
    "ax[0].set_xlabel(r\"Time $t$\",fontsize=18)\n",
    "ax[0].set_ylabel(r\"$x$ and $y$\",fontsize=18)\n",
    "ax[0].legend(loc=\"best\",fontsize=12)\n",
    "ax[0].set_title(\"Cells that choose fate 1\",fontsize=20)\n",
    "\n",
    "ax[1].plot(TimeRange,avg_traj[2,0],c=CB_colors[1],label=\"x\")\n",
    "ax[1].fill_between(TimeRange,avg_traj[2,0]-std_traj[2,0],avg_traj[2,0]+std_traj[2,0],color=CB_colors[1],alpha=0.5)\n",
    "ax[1].plot(TimeRange,avg_traj[2,1],c=CB_colors[2],label=\"y\")\n",
    "ax[1].fill_between(TimeRange,avg_traj[2,1]-std_traj[2,1],avg_traj[2,1]+std_traj[2,1],color=CB_colors[2],alpha=0.5)\n",
    "ax[1].axhline(thresh[0],color='k',ls='--')\n",
    "ax[1].set_xlabel(r\"Time $t$\",fontsize=18)\n",
    "ax[1].legend(loc=\"best\",fontsize=12)\n",
    "ax[1].set_title(\"Cells that choose fate 2\",fontsize=20)\n",
    "\n",
    "fig.suptitle(fr\"Average trajectories\",fontsize=22,y=1.02)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aad1bf32",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now let's vary the constant external bias, keeping theta and noise fixed:\n",
    "ext_bias_arr = np.array([i/100 for i in range(11)])/T\n",
    "fate_frax_bias = np.zeros((len(ext_bias_arr),4))\n",
    "\n",
    "for ext_bias_idx in range(len(ext_bias_arr)):\n",
    "#     print(f\"field: {ext_bias_arr[ext_bias_idx]}\")\n",
    "#     print(f\"time integral of field: {np.trapz(ext_bias_arr[ext_bias_idx]*np.ones_like(TimeRange),x=TimeRange)}\\n\")\n",
    "    \n",
    "    #Generating the constant external field:\n",
    "    field = ext_bias_arr[ext_bias_idx]\n",
    "    \n",
    "    #Initial conditions, exactly at (0,0):\n",
    "    x0 = np.zeros((2,num_traj))\n",
    "    \n",
    "    #Time evolution given the constant external field:\n",
    "    traj_ext_bias = np.empty((2,num_traj,Nt))\n",
    "    traj_ext_bias = Solver_const_ext_bias(x0,a,r,k,n,theta,alpha,field)\n",
    "    \n",
    "    #Calculating stabilities and fate fractions:\n",
    "    stability_factrs,fates = fate_classifier(traj_ext_bias)\n",
    "    fate_frax_bias[ext_bias_idx] = fate_fractions(fates)\n",
    "#     avg_traj,std_traj = traj_moments(traj_ext_bias,fates)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f876ea97",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(ext_bias_arr,fate_frax_bias[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc0dc9f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fate_frax_bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e5bf3a7d",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'X' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [6]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mX\u001b[49m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'X' is not defined"
     ]
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "029c3feb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}