coral_3D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c145ef165f1d4877a9cf4ad9227d84bc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=9, description='i', max=19), Output()), _dom_classes=('widget-interact',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "from ipywidgets import interact\n",
    "\n",
    "# Size of the 3D grid\n",
    "N = 20\n",
    "\n",
    "# Time step\n",
    "dt = 0.05\n",
    "\n",
    "# Initialize U and V\n",
    "U = np.full((N, N, N), 0.0)  # Initial condition for U\n",
    "V = np.zeros((N, N, N))  # Initial condition for V\n",
    "\n",
    "# Initialize small area in center of V = 1\n",
    "center_size = 2\n",
    "V[N // 2 - center_size:N // 2 + center_size,\n",
    "  N // 2 - center_size:N // 2 + center_size,\n",
    "  N // 2 - center_size:N // 2 + center_size] = 1\n",
    "\n",
    "# Parameters\n",
    "Du, Dv = 1.0, 0.5  # Adjusted parameters for the coral pattern\n",
    "\n",
    "# Initialize F and k\n",
    "F_init = 0.0545\n",
    "k_init = 0.062\n",
    "\n",
    "# Function to compute F\n",
    "def compute_F(F_init, x, y, z, t):\n",
    "    # Compute the square distance from the center of the grid\n",
    "    d_squared = (x - N/2)**2 + (y - N/2)**2 + (z - N/2)**2\n",
    "\n",
    "    # Initial Gaussian-like function for the spatial distribution\n",
    "    spatial_dependence = np.exp(-d_squared / (2.0 * (N/4)**2))\n",
    "\n",
    "    # Time dependence: the concentration decreases at the center and increases in outer locations over time\n",
    "    time_dependence = 1 + np.tanh((d_squared - (N*t/2000)**2)/(N**2/16))\n",
    "\n",
    "    # Combine the spatial and temporal dependence\n",
    "    F = F_init * spatial_dependence * time_dependence\n",
    "\n",
    "    return F\n",
    "\n",
    "\n",
    "# Function to compute k\n",
    "def compute_k(k_init, x, y, z, t):\n",
    "    # k is modulated by the distance from the origin (increasing with time), \n",
    "    # creating a moving wave pattern\n",
    "    k = k_init + 0.0100 * np.sin(np.sqrt((x - t)**2 + (y - t)**2 + (z - t)**2))\n",
    "    return k\n",
    "\n",
    "# Function to calculate Laplacian using convolution\n",
    "def laplacian(Z):\n",
    "    laplacian_kernel = np.array([\n",
    "        [[0.05, 0.2, 0.05],\n",
    "         [0.2, -1, 0.2],\n",
    "         [0.05, 0.2, 0.05]],\n",
    "\n",
    "        [[0.05, 0.2, 0.05],\n",
    "         [0.2, -1, 0.2],\n",
    "         [0.05, 0.2, 0.05]],\n",
    "        \n",
    "        [[0.05, 0.2, 0.05],\n",
    "         [0.2, -1, 0.2],\n",
    "         [0.05, 0.2, 0.05]]\n",
    "    ])\n",
    "\n",
    "    return signal.convolve(Z, laplacian_kernel, mode='same')\n",
    "\n",
    "# Store each frame of the simulation\n",
    "frames_V = []\n",
    "frames_F = []\n",
    "frames_k = []\n",
    "\n",
    "# Run the simulation\n",
    "for i in range(2000):  # 2000 steps for computational efficiency\n",
    "    # Compute F and k\n",
    "    x, y, z = np.indices(U.shape)\n",
    "    t = np.full(U.shape, i)  # Expand i to match the shape of indices\n",
    "\n",
    "    F = compute_F(F_init, x, y, z, t)\n",
    "    k = compute_k(k_init, x, y, z, t)\n",
    "\n",
    "    # Update U and V with the constraint on surface area\n",
    "    U_new = U + dt * (Du * laplacian(U) - U*V*V + F*(1-U))\n",
    "    V_new = V + dt * (Dv * laplacian(V) + U*V*V - (F+k)*V)\n",
    "\n",
    "    # Clip values to keep within the range of 0 and 1\n",
    "    U_new = np.clip(U_new, 0, 1)\n",
    "    V_new = np.clip(V_new, 0, 1)\n",
    "\n",
    "    # Update U and V\n",
    "    U, V = U_new, V_new\n",
    "\n",
    "    if i % 100 == 0:  # Store every 100th frame\n",
    "        frames_V.append(V.copy())\n",
    "        frames_F.append(F.copy())\n",
    "        frames_k.append(k.copy())\n",
    "\n",
    "# Function to display a frame\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
    "plt.close()\n",
    "\n",
    "def display_frame(i):\n",
    "    fig.clear()\n",
    "    \n",
    "    ax1 = fig.add_subplot(131, projection='3d')\n",
    "    ax1.voxels(frames_V[i] > 0.2, facecolors='b', alpha=0.2)\n",
    "    ax1.set_title('Frame: %d' % i)\n",
    "\n",
    "    ax2 = fig.add_subplot(132, projection='3d')\n",
    "    ax2.scatter(*np.indices(frames_F[i].shape), c=frames_F[i].ravel(), alpha=0.5)\n",
    "    ax2.set_title('F Distribution')\n",
    "\n",
    "    ax3 = fig.add_subplot(133, projection='3d')\n",
    "    ax3.scatter(*np.indices(frames_k[i].shape), c=frames_k[i].ravel(), alpha=0.5)\n",
    "    ax3.set_title('k Distribution')\n",
    "\n",
    "    display(fig)\n",
    "\n",
    "# Use a slider to navigate through the frames\n",
    "interactive_plot = interact(display_frame, i=(0, len(frames_V)-1))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rd-corals",
   "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.6"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}