{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction to Agent-Based Modeling (Python-based)\n", "## Paul Macklin, Indiana University (macklinp@iu.edu) \n", "## Version 2: February 2, 2022\n", "\n", "This Jupyter notebook accompanies a two-lecture series that introduces agent-based modeling (in biology) by building a basic agent-based model in Python, and then explores more developed cloud-hosted agent-based models created with [PhysiCell](http://PhysiCell.org). \n", "\n", "This notebook and the accompanying lectures are available as open education resources under the [CC-BY 4.0 license](https://creativecommons.org/licenses/by/4.0/). Please visit the following GitHub repository for the latest version: \n", "\n", "https://github.com/physicell-training/Python_ABM_crashcourse\n", "\n", "### Version History\n", "#### Version 1. February 16, 2021\n", "Initial version. \n", "\n", "#### Version 2. February 2, 2022\n", "Fixed typographical errors, added improved plotting, and fixed nuisance np.int vs int warnings. \n", "\n", "### Section 1. Declaring a class in Python \n", "Let's start by declaring a class. \n", "\n", "The name of this class is `Agent`. It is a template for a software object. Any specific software object created from this template is an *instance* of the class. \n", "\n", "Its default constructor is `__init__`. Its single input variable (argument) is **self**: the software object (instance) created by the class. Inside of `__init__`, we define the member data, which are attached to `self`. \n", "\n", "Every Python class must have a default constructor with syntax `__init__(self)`. \n", "\n", "This particular class has a member function (`move_me`), with arguments `self` (always required!), `dx`, and `dy`. It also has a method `display` to show its current state. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "class Agent:\n", " def __init__( self ): # default constructor\n", " self.hidden_variable = False; # no such thing as private in Python\n", " self.position = [0,0]; \n", " self.hunger = 0.0;\n", " self.ID = 0; \n", " return; \n", "\n", " def move_me( self , dx, dy ): # member function \n", " self.position[0] += dx ; \n", " self.position[1] += dy; \n", " return; \n", " \n", " def display( self ):\n", " print( str(self.ID) + ' at ' + str(self.position) + ' has hunger level ' \n", " + str(self.hunger) )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's create an instance of this class named Bob." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 at [0, 0] has hunger level 0.0\n" ] } ], "source": [ "Bob = Agent(); \n", "Bob.display() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's set Bob's ID to 1. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 at [0, 0] has hunger level 0.0\n" ] } ], "source": [ "Bob.ID = 1 \n", "Bob.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's use the member function to move Bob one right, one down: " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 at [2, -2] has hunger level 0.0\n", "1 at [2, -2] has hunger level 0.0\n", "2 at [0, 0] has hunger level 0.5\n" ] } ], "source": [ "Bob.move_me( 1 , -1 )\n", "Bob.display()\n", "\n", "Alice = Agent()\n", "Alice.ID = 2\n", "Alice.hunger = 0.5 \n", "\n", "Bob.display()\n", "Alice.display()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 2. Loading key libraries. \n", "Let's load the \"standard\" Python libraries we'll need, and set matplotlib to plot inline (for now). " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt \n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 3. Creating a basic microenvironment \n", "\n", "Let's now write a basic `Environment` class.\n", "\n", "We'll introduce a basic chemical factor (a growth substrate) $c$ that obeys diffusion and decay: \n", "\n", "$$ \\frac{ \\partial c}{\\partial t} = D \\nabla^2 c - \\lambda c $$ \n", "\n", "First, we're going to add a basic Cartesian mesh: $X$ will hold regularly-spaced x-coordinates $a \\le x \\le b$, and $Y$ will hold y-coordinates with $c \\le y \\le d$. We'll have $m \\times n$ lattice points, with \n", "\n", "$$ \\Delta x = \\frac{b-a}{m-1} , \\textrm{ and } \\Delta y = \\frac{d-c}{n-1}.$$\n", "\n", "And then $x_i = a + i\\Delta x$ and $y_j = c+j \\Delta y$. \n", "\n", "We'll use numpy linspaces for X and Y, and a numpy matrix (initially zeros) to store the substrates. \n", "\n", "To make this easy, we'll put most of the initial sizing into the default constructor (`__init__`), and we'll use the `setup` function to set the physical parameters and apply the initial condition. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Declare environment class\n", "class Environment:\n", " def __init__( self , shape=[-100,-100,100,100], m=2 , n=2 ):\n", " self.a = shape[0]\n", " self.b = shape[2]\n", " self.c = shape[1]\n", " self.d = shape[3]\n", " self.m = m; \n", " self.n = n; \n", " self.dx = (self.b-self.a)/(m-1); \n", " self.dy = (self.d-self.c)/(n-1);\n", " self.X = np.linspace( self.a, self.b , m )\n", " self.Y = np.linspace( self.c, self.d, n )\n", " self.C = np.zeros((m,n)); \n", " return; \n", " def setup( self , D=1000, decay=0.1 , initial=1.0, boundary=1.0 ):\n", " self.D = D; \n", " self.decay = decay; # In Python you can evidently declare more class elements in methods\n", " # set boundary values \n", " self.C = boundary * np.ones((self.m,self.n))\n", " # set interior values to initial value \n", " for j in range(1,self.n-1):\n", " for i in range(1,self.m-1):\n", " self.C[i,j] = initial ; \n", " print( 'Length scale: ' + str( np.sqrt(self.D/self.decay) ) )\n", " return; " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's test it on a $[-100,200] \\times [0,400]$ domain with $31\\times 41$ lattice points. Diffusion coefficient 1000, decay rate 0.1, initial value of 0.5, and boundary value 1. \n", "\n", "We'll do a contour plot to make sure it's good. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length scale: 100.0\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD8CAYAAACRvtrKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVp0lEQVR4nO3de7BdZX3G8e9jSISUAEmBGEEhMiESHRMxIpfWUhCN8RLoSIUOmGFw0BYQnVoNOq12Os5gBQVFYCIgQRAaI0rKREhMRbTKnQCBQAmExoSYcJVLuSX8+sd695uVk31Z5+xzctYOz2dmz977XWvttdbss35n3fb7KCIwMwN4w3AvgJnVhwuCmWUuCGaWuSCYWeaCYGaZC4KZZZULgqQRku6SdF16P07SEkkPpeexpXHPlLRS0oOSPjQUC272eiDpUkkbJC1vMVySvpu2t3skHdhn+BbbbSf92UM4A1hRej8HWBoRk4Cl6T2SpgDHAe8AZgAXSBrRj/mY2WaXUWxHrXwYmJQepwAX9hned7ttq1JBkLQ38BHg4lLzLGBeej0POLrUfnVEvBwRq4CVwEFVF8jMNouIm4Cn2owyC7g8CjcDu0maAC2327Z2qDjeucCXgDGltvERsS4t9DpJe6b2vYCbS+OtSW1bkHQKRUVjp9F6z8T9qi6K2dC5/95Xn4iIPbr5jMMO3zGeeeq1qvO7D3ip1DQ3Iub2Y3Z7AX8ovW9sb+tovt221XErlPRRYENE3CHp8AqfqSZtW90fnVZ6LsA73jUqfnzd+K0m+tHTh1aYndnAnDj2d1u1Tdtnzf92+7nPPPUazf6em5m2z5qXImJ6F7Nrur0NYLsFqu0hHAZ8XNJMYEdgF0lXAOslTUh7BxOADWn8NcBbStPvDTxWdYEafvT0oSxePbm/k5n1S7Oi0GNabW+foMl2GxEntPuwjgUhIs4EzgRIleaLEXGCpG8Bs4Gz0vO1aZKFwI8lfRt4M8XJjlurrh1sLgbPr9q1P5OZ9ctiin84PV4UFgKnSboaeB/wp3Qo33S77fRh3Ry4nwXMl3QysBo4FiAi7pM0H7gf2AicGhGbqn5ouRiMedi3SdjQeY5da18UJF0FHA7sLmkN8DVgJEBEXAQsAmZSnLz/P+CkbubXr4IQETcCN6bXTwJHthjvG8A3+rsw5WLw1kUb+zu5WT+Nqn1RiIjjOwwP4NQO49xI2m47qc2p/WZ7BjstXzvMS2Xbt70oF4U6FoRtzfvkZpa5IJhZ5oJgZpkLgpllLghmlrkgmFnmgmBmmQuCmWUuCGaWuSCYWeaCYGaZC4KZZS4IZpa5IJhZ5oJgZpkLgpllLghmlnUsCJJ2lHSrpLsl3SfpX1P71yWtlbQsPWaWpnGUm1kPqtKF2svAERHxvKSRwG8l/SIN+05EnF0euU+U25uBX0ravz8drZrZ8Oi4h5Aiop5Pb0emx1bBKyWOcjPrUVWzHUdIWkYRxrIkIm5Jg05LibOXltKfW0VLmVnNVSoIEbEpIqZRpMIcJOmdFCmz+wHTKHLkzkmjV4pyk3SKpNsl3f50xRw8s9cjSTPS+biVkuY0GT5W0s/SP+db0/bZGLabpAWSHpC0QtIh7ebVr6sMEfEMRf/uMyJifSoUrwE/YPNhQaUot4iYGxHTI2L62HG+2GHWjKQRwPcpYt+nAMen83RlXwGWRcS7gE8B55WGnQdcHxFvB6bSIRq+ylWGPSTtll7vBHwAeKAROZ0cAyxPrxcCx0l6o6SJDCDKzcyyg4CVEfFIRLwCXE1xnq5sCrAUICIeAPaVNF7SLsD7gUvSsFfSP/WWqlxlmADMS5XqDcD8iLhO0o8kTaM4HHgU+EyaaVdRbma97slNO/cjuXz+7pJuLzX0jYNvdk7ufX0+5G7gbyiuAB4E7EOxZ74JeBz4oaSpwB3AGRHxQqulqRL2eg/w7ibtJ7aZZkBRbmavQ090iIOvck7uLOC8dOL/XuAuin/GI4EDgdMj4hZJ5wFzgH9uNbPaRLmZWVMdz8lFxLOkkFdJAlalx2hgTemq4AKKgtCSz+aZ1dttwCRJEyWNorjpb2F5hHQlYVR6+2ngpoh4NiL+CPxB0uQ07EiKQ/mWvIdgVmMRsVHSacANwAjg0nSe7rNp+EXAAcDlkjZRbPAnlz7idODKVDAeoUNcvAuCWc1FxCJgUZ+2i0qvf09xNa/ZtMuAducotuBDBjPLXBDMLHNBMLPMBcHMMhcEM8tcEMwsc0Ews8wFwcwyFwQzy1wQzCxzQTCzzAXBzDIXBDPLXBDMLHNBMLOsm2zHcZKWSHooPY8tTeNsR7MeVGUPoZHtOJUilGWGpIMp+mZbGhGTKLqAngNbZTvOAC5IPTabWc11k+04C5iX2ucBR6fXznY061HdZDuOj4h1AOl5zzR6pWxHR7mZ1U832Y6tVMp2dJSbWf0MONsRWN+Ic0vPG9JolbIdzax+BpztSNE3/Ow02mzg2vTa2Y5mParKHsIE4FeS7qEIjVgSEddRxEcdJekh4Kj0noi4D2hkO16Psx3NulIhDn5XSf9ZujXgpNKwL6S25ZKukrRju3l1k+34JEUSTLNpnO1oNghKcfBHURyO3yZpYUSUE5hOBe6PiI9J2gN4UNKVwB7A54ApEfFiCmE+Dris1fx8Ns+s3qrEwQcwJuU67gw8RRH2CsU//Z0k7UCR9dj2fJ6Tm8wG2bOv7Mji1ZM7j1gYjDj48ynO3T0GjAE+GRGvAWslnQ2sBl4EFkfE4nYL44JgNrwGIw7+Q8Ay4AhgP2CJpN9QZEHOAiYCzwA/kXRCRFzRamY+ZDCrtyqX8U8Crkl3Fa+kiIJ/O8UVwVUR8XhEvApcAxzabmYuCGb11jEOnuKQ4EgASeOByRRJz6uBgyWNTucXjgRWtJuZDxnMaqxiHPy/AZdJupfiEOPLEfEE8ISkBcCdFCcZ7wLmNptPgwuCWc1ViIN/DPhgi2m/Bnyt6rx8yGBmmQuCmWUuCGaWuSCYWeaCYGaZC4KZZS4IZpa5IJhZ5oJgZpkLgpllLghmllXpZPUtkn4laUXqm+2M1P51SWslLUuPmaVpHOVm1oOq/LhpI/CPEXGnpDHAHZKWpGHfiYizyyP3iXJ7M/BLSfu7o1Wz+qsS5bYuIu5Mr5+j+D31VklMJY5yM+tR/TqHIGlfih6Yb0lNp0m6R9KlpfRnR7mZ9ajKBUHSzsBPgc9HxLPAhRT9t00D1gHnNEZtMrmj3Mx6QNWw15EUxeDKiLgGICLWp8zH14AfsPmwwFFuZj2qylUGAZcAKyLi26X2CaXRjgGWp9eOcjPrUVWuMhwGnAjcmyLhAb4CHC9pGsXhwKPAZ6CIcksJMfdTXKFwlJtZj6gS5fZbmp8XWNSkrTGNo9zMepDP5plZ5oJgZpkLglnNdRkH33bavlwQzGqsFAf/YWAKxcn8KX1Ga8TBTwUOB86RNKritFtwQTCrt27i4KtMuwUXBLN6q/JTgPOBAyhuALwXOCPdMFjpZwRljnIzG2SbXh7B86t2rTr67pJuL72fGxHl/MVu4uAr/YygzAXBbHg9ERHT2wyvGgd/VkQEsFJSIw6+3z8j8CGDWb11EwdfZdoteA/BrMa6jIOn2bTt5ueCYFZzXcbBbzVtOz5kMLPMBcHMMhcEM8tcEMwsc0Ews8wFwcwyFwQzy7qJchsnaYmkh9Lz2NI0jnIz60FV9hAaUW4HAAcDp6bfVM8BlkbEJGBpet83ym0GcEH6XbaZ1Vw3UW6zgHlptHnA0em1o9zMelQ3UW7jI2IdFEUD2DON5ig3sx7VTZRby1GbtDnKzawHDDjKDVjfSG9KzxtSu6PczHrUgKPcKH5XPTu9ng1cW2p3lJtZD+omyu0sYL6kkyk6aDgWHOVm1su6iXKD1EtLk2kc5WbWg3w2z8wyFwQzy1wQzCxzQTCzzAXBzDIXBDPLXBDMaq5CHPw/SVqWHsslbUrdEzTtuqAdFwSzGqsS6R4R34qIaRExDTgT+HVENBKgm3Vd0JILglm99TfS/XjgKmjbdUFLLghm9VY50l3SaIpOiX7aZNi+bO66oCVHuZkNshEvw5iHK/+vHYw4+IaPAf+dDhc2f0D1rgtcEMyG2WDEwTccRzpcaGjRdUFLPmQwq7dKke6SdgX+is3dELTruqAlFwSzGouIjUAj0n0FML8RB9+IhE+OARZHxAultkbXBUeULkvObDc/HzKY1VynOPj0/jLgsj5t7bouaMp7CGaWuSCYWeaCYGZZlU5WL5W0QdLyUtvXJa1tdqLCMW5mvavKHsJlFHc/9fWdxv3T6aSHY9zMelyVKLebgKc6jZc4xs2sh3VzDuE0SfekQ4pG8nN/7rt2lJtZzQy0IFwI7AdMA9YB56T2yvddO8rNrH4GtCVGxPqI2BQRrwE/YPNhgWPczHrYgApCI9MxOQZoXIFwjJtZD+t467Kkq4DDKX6muQb4GnC4pGkUhwOPAp8Bx7iZ9boqUW7HN2m+pM34jnEz61E+m2dmmQuCmWUuCGaWuSCYWeaCYGaZC4KZZS4IZpa5IJhZ5oJgZpkLgpllLghmlrkgmNWcpBmpj9KVkua0GOfw1L/pfZJ+3WfYCEl3Sbqu07wc1GJWY6lP0u8DR1H0N3KbpIURcX9pnN2AC4AZEbFa0p59PuYMitSnXTrNz3sIZvV2ELAyIh6JiFeAqyn6Li37O+CaiFgNEBEbGgMk7Q18BLi4ysy8h2A2yEa8FIx98JWqo3eKg2/WT+n7+nzG/sBISTcCY4DzIuLyNOxc4EupvSMXBLPh1SkOvko/pTsA7wGOBHYCfi/pZopCsSEi7pB0eJWFcUEwq7cq/ZSuoSgsLwAvSLoJmAocCHw8BSntCOwi6YqIOKHVzHwOwazebgMmSZooaRRFENLCPuNcC/ylpB0kjaY4pFgREWdGxN4RsW+a7r/aFQMYeJTbOElLJD2UnseWhjnKzWyQRMRG4DTgBoorBfNT36WflfTZNM4K4HrgHopOjS+OiOWtPrOdgUa5zQGWRsQkYGl67yg3syEQEYsiYv+I2C/1WUpEXBQRF5XG+VZETImId0bEuU0+48aI+GineQ00ym0WMC+9ngccXWp3lJtZjxroOYTxEbEOID03boRwlJtZDxvsk4qOcjPrYQPdEtc30pvSc+POKEe5mfWwgRaEhcDs9Ho2xWWPRruj3Mx61ECj3M4C5ks6GVgNHAuOcjPrdQONcoPiNslm4zvKzaxH+WyemWUuCGaWuSCYWeaCYGaZC4KZZS4IZpa5IJhZ5oJgZpkLgpllLghmlrkgmFnmgmBmmQuCmWUuCGaWuSCYWeaCYFZzneLgUxT8n1Ic/DJJ/1IatpukBZIekLRC0iHt5uUoN7MaqxIHn/ymRe7CecD1EfGJlPw0ut38vIdgVm9V4uCbkrQL8H7gEoCIeCUinmk3TVd7CJIeBZ4DNgEbI2K6pHHAfwD7Ao8CfxsRT3czH7Ne8oYXX2Wn5Wurjj4YcfAAh0i6m6KX8y9GxH3A24DHgR9KmgrcAZyRQmGbL3vVpW7jryNiWinSumnMm5k19UQjnyQ95vYZXiXr5E5gn4iYCnwP+Hlq34EiAfrCiHg38AIdtsehOGRoFfNmZv3XMeskIp6NiOfT60XASEm7p2nXRMQtadQFFAWipW4LQgCLJd0h6ZTU1irmzcz6r2McvKQ3SVJ6fRDFdv1kRPwR+IOkyWnUIykiElrq9irDYRHxmKQ9gSWSHqg6YSogpwBM2MsB0WbNRMRGSY04+BHApY04+DT8IuATwN9L2gi8CBwXEY3DitOBK1MxeQQ4qd38uioIEfFYet4g6WcUZ0TXS5oQEev6xLz1nXYuMBfgHe8a1TT/0czyYcCiPm3lKPjzgfNbTLsMmN5sWDMDPmSQ9GeSxjReAx8EltM65s3Maq6bPYTxwM/SocsOwI8j4npJt9Ek5s3M6m/ABSEiHgGmNml/khYxb2ZWb75T0cwyFwQzy1wQzCxzQTCzzAXBzDIXBDPLXBDMLHNBMLPMBcHMMhcEM8tcEMwsc0Ews8wFwcwyFwQzy1wQzCxzQTCzrDZRbieO/R0Ai5kMD48d5qWx14OnJ4/iuf1eY+eJf+KDb31wuBenFmpTEKBUFD4wmacZy9OTJw7zEtn2rFwMGn97r3e1KgiwZVF4ftWuw7w0tj1zMdjakBUESTMokmdHABdHxFlVp93i8MFsiPRKMai6LUl6L3Az8MmIWJDavgB8miJU6V7gpIh4qdW8hqQg9CPCuqVe+KKst/XC31jVbSmN902KQJdG217A54ApEfGipPkUyU+XtZrfUO0h5AjrtGCNCOvKBQF64wszG2JVt6XTgZ8C7+3TvgOwk6RXgdH0yYXsa6gKQscI63KUG/DytH3WLB+iZRluuwNPDPdCDIHtdb26Pk599tUNN1y/9nu7Vxx9x27j4NOewDHAEZQKQkSslXQ2RT7Ki8DiiFjcbmGGqiB0jLAuR7lJur0UJ79d2V7XbXter24/IyJmDMayJFXi4M8FvhwRm1JwUjGhNJZib2Ii8AzwE0knRMQVrWY2VAWhY4S1mVVSZVuaDlydisHuwMwU/DoSWBURjwNIugY4FGhZEIbqTsWOEdZmVknHbSkiJkbEvhGxL7AA+IeI+DnFocLBkkanuPgjgRXtZjYkewitIqzbTDK3zbBet72um9drG6gYB99q2lskLQDuBDYCd9Fh/bQ5Rt7MXu/84yYzy1wQzCzb5gVB0rGS7pP0mqTpfYadKWmlpAclfajU/h5J96Zh31X52kpNSZqR1mOlpDnDvTz9JelSSRskLS+1jZO0RNJD6XlsaVjT765uJL1F0q8krUh/h2ek9p5ft0EREdv0ARxAcfPHjcD0UvsU4G7gjRTXTR8GRqRhtwKHUFyT/QXw4W293P1cxxFp+d8GjErrNWW4l6uf6/B+4EBgeant34E56fUc4Judvru6PYAJwIHp9Rjgf9Ly9/y6DcZjm+8hRMSKiGj24/NZwNUR8XJErAJWAgdJmgDsEhG/j+Ibuhw4etst8YDk200j4hWgcbtpz4iIm4Cn+jTPAual1/PY/D00/e62xXL2V0Ssi4g70+vnKC7D7cV2sG6DoU7nEJrdorlXeqxp0l5nrdal142PiHVQbFjAnqm9J9dX0r7Au4Fb2M7WbaCG6teOvwTe1GTQVyPi2laTNWmLNu111ovL3I2eW19JO1P8GOjzEfFsm9NSPbdu3RiqG5M+MIDJWt2iuSa97tteZ9vrrdvrJU2IiHXpUG5Dau+p9ZU0kqIYXBkR16Tm7WLdulWnQ4aFwHGS3ihpIjAJuDXtvj0n6eB0deFTQKu9jLrYXm/dXgjMTq9ns/l7aPrdDcPydZT+hi4BVkTEt0uDen7dBsUwnOU9hqLqvgysB24oDfsqxVncByldSaD48cbyNOx80h2WdX4AMynOYD9Mcag07MvUz+W/ClgHvJq+r5OBPweWAg+l53Gdvru6PYC/oNjlvwdYlh4zt4d1G4yHb102s6xOhwxmNsxcEMwsc0Ews8wFwcwyFwQzy1wQzCxzQTCz7P8BeLla2B0TmAMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "E = Environment( [-100,0,200,400] , 31,41 ) \n", "\n", "E.setup( 1000, 0.1 , 0.5 , 1 )\n", "\n", "plt.contourf( E.X, E.Y, np.transpose( E.C ) )\n", "plt.axis('image')\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 4. Adding diffusion and built-in plotting \n", "We'll solve the diffusion-decay PDE with an explicit finite difference method. If $(x_i,y_j)$ \n", "is some lattice site and $t_n = n \\Delta t$ is a time, then let's write: \n", "$$\n", "c^n_{i,j} = c( x_i, y_j , t_n )\n", "$$\n", "\n", "Using this notation and standard finite differenes, we have: \n", "$$\n", "\\frac{ c^{n+1}_{i,j} - c^{n}_{i,j} }{\\Delta t} = \n", "D \\left( \\frac{ c^n_{i+1,j}-2 c^n_{i,j} - c^n_{i-1,j} }{\\Delta x^2 } + \n", "\\frac{ c^n_{i,j+1}-2 c^n_{i,j} - c^n_{i,j-1} }{\\Delta y^2 } \n", "\\right) - \\lambda c^n_{i,j}\n", "$$\n", "\n", "Then, we can algebraically solve for $c^{n+1}_{i,j}$: \n", "$$\n", "c^{n+1}{i,j} = \n", "c^{n}{i,j} \n", "+ \n", "\\frac{D \\Delta t}{\\Delta x^2} \n", "\\left( c^n_{i+1,j}-2 c^n_{i,j} - c^n_{i-1,j} \\right)\n", "+ \n", "\\frac{D \\Delta t}{\\Delta y^2} \n", "\\left( c^n_{i,j+1}-2 c^n_{i,j} - c^n_{i,j-1} \\right) \n", "- \\lambda \\Delta t c^{n}_{i,j}\n", "$$\n", "We can clean this up to: \n", "$$\n", "A = \\frac{D\\Delta t}{\\Delta x^2} , \\hspace{0.25in} \n", "B = \\frac{D \\Delta t}{\\Delta y^2} , \\hspace{0.25in}\n", "C = \\Delta t \\lambda \n", "$$\n", "and then \n", "$$\n", "c^{n+1}_{i,j} = \\left( 1 - (2A+2B) - C \\right)c^n_{i,j} \n", "+ A c^n_{i-1,j} + A c^n_{i+1,j} + B c^n_{i,j-1} + B c^n_{i,j+1} \n", "$$\n", "We'll iterate through our interior (non-boundary) points using this computational stencil. \n", "\n", "One more thing: when you do `b = a` in Python, it doesn't actually copy `a` to create a new `b`. It instead makes a pointer to `a` and stores it at `b`. This makes it efficient to pass objects to functions, but it has a serious downside: operating on your \"copy\" `b` changes the original `a`. You have to use `copy` (or sometimes `deepcopy`) to get around this Python design choice. \n", "\n", "Let's re-define the `Environment` class with this algorithm. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import copy \n", "\n", "# Declare environment class\n", "class Environment:\n", " def __init__( self , shape=[-100,-100,100,100], m=3 , n=3 ):\n", " self.a = shape[0]\n", " self.b = shape[2]\n", " self.c = shape[1]\n", " self.d = shape[3]\n", " self.m = m; \n", " self.n = n; \n", " self.dx = (self.b-self.a)/(m-1); \n", " self.dy = (self.d-self.c)/(n-1);\n", " self.X = np.linspace( self.a, self.b , m )\n", " self.Y = np.linspace( self.c, self.d, n )\n", " self.C = np.zeros((m,n)); \n", " return; \n", " def setup( self , D=1000, decay=0.1 , initial=1.0, boundary=1.0 ):\n", " self.D = D; \n", " self.decay = decay; \n", " # set all to the boundary values\n", " self.C = boundary*np.ones((self.m,self.n))\n", " # now set interior values to the initial value \n", " for j in range(1,self.n-1):\n", " for i in range(1,self.m-1):\n", " self.C[i,j] = initial ; \n", " print( 'Length scale: ' + str( np.sqrt(self.D/self.decay) ) )\n", " return; \n", " def update( self , dt=0.001 ):\n", " # define constants for simplicity \n", " A = self.D * dt / ( self.dx**2 )\n", " B = self.D * dt / ( self.dy**2 )\n", " C = self.decay * dt \n", " # copy the prior solution. \n", " old = self.C.copy(); # make sure this is a real copy and not a reference. \n", " # finite differences for decay-diffusion \n", " for i in range(1,self.m-1):\n", " for j in range(1,self.n-1):\n", " self.C[i,j] = (1-2*A-2*B-C)*old[i,j] + A*(old[i-1,j]+old[i+1,j]) + B*( old[i,j-1]+old[i,j+1] ); \n", " return; \n", " def plot( self ):\n", " plt.clf()\n", " plt.contourf( self.X, self.Y, np.transpose( self.C ) )\n", " plt.axis('image')\n", " plt.colorbar()\n", " plt.xlabel('x')\n", " plt.ylabel('y')\n", " plt.show()\n", " plt.pause(0.0001) # helps with animation\n", " \n", "# a startup function for full-screen plots \n", "def fullscreen_setup( fignum=1 , pause_time=2):\n", " %matplotlib qt \n", " plt.figure(fignum)\n", " plt.clf()\n", " figManager = plt.get_current_fig_manager()\n", " figManager.full_screen_toggle() \n", " plt.pause(pause_time); # some time to switch windows if needed " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's try it out. We'll define a domain, with initial condition 0.5 and boundary value 1. \n", "\n", "Let's use a mesh size of $\\Delta x = \\Delta y = 20.$ For a finite difference scheme like this, there are restrictions on how big of a time step size we can use without losing accuracy and numerical stability: \n", "\n", "$$\\Delta t \\le \\frac{ \\Delta x^2}{2n D},$$\n", "\n", "where $n$ is the number of spatial dimensions (here, $n = 2$). \n", "\n", "If $\\Delta x = 20$ and $D = 10000$, then we need: \n", "$$\\Delta t \\le \\frac{400}{4\\cdot 10000} = 0.01$$\n", "\n", "Let's run update every 0.01 minutes and plot every 200 time steps. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length scale: 1000.0\n", "0.01\n" ] } ], "source": [ "# declare an environment \n", "E = Environment( [-250,-250,250,250] , 26,26 ) \n", "# set up\n", "E.setup( 10000, 1 , 0.5 , 1 )\n", "\n", "dt = E.dx**2 / (2 * 2 * E.D )\n", "\n", "print(dt)\n", "\n", "# prepare for full-screen plots\n", "fullscreen_setup()\n", "\n", "# main test loop\n", "for n in range(100):\n", " if( n % 10 == 0 ):\n", " E.plot()\n", " # plt.show() \n", " # plt.pause(0.001)\n", " E.update( dt )\n", " \n", "E.plot() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 5. A cell class (version 1)\n", "Let's write a cell class that does the following: \n", "* updates cell velocity by exchanging spring-like forces with other cells. \n", "* there is a random component to velocity (due to random actin dynamics etc)\n", "* updates cell position based on the velocity. \n", "We'll also make a list of cells to keep track of them all. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "all_cells = list(); \n", "\n", "class Cell:\n", " def __init__( self ):\n", " self.position = np.array( [0.0,0.0] ); \n", " self.velocity = np.array( [0.0,0.0] );\n", " self.mechanics_distance = 30.0; \n", " self.equilibrium_spacing = 20.0;\n", " self.mechanics_contant = 1.0; \n", " def update_velocity( self, dt , env , all_cells ):\n", " self.velocity = np.array([0.0,0.0])\n", " # spring-like interactions with cells \n", " for c in all_cells: \n", " displacement = c.position - self.position \n", " dist = np.linalg.norm( displacement )\n", " if( dist < self.mechanics_distance and dist > 1e-16 ):\n", " dv = dist - self.equilibrium_spacing \n", " displacement = displacement / dist; \n", " displacement = displacement * dv; \n", " displacement = displacement * self.mechanics_contant; \n", " self.velocity = self.velocity + displacement \n", " angle = 6.28318530718 * np.random.uniform(); \n", " perturbation_size = 0.1 * self.mechanics_contant; \n", " self.velocity[0] += perturbation_size * np.cos(angle)\n", " self.velocity[1] += perturbation_size * np.sin(angle)\n", " return\n", " def update_position( self, dt , env, all_cells ):\n", " self.position = self.position + dt*self.velocity; \n", " return; " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's make a few cells with random positions contained in the environment and plot. " ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "all_cells = list()\n", "\n", "# make a plotting function \n", "def plot_cells( env, all_cells ):\n", " env.plot()\n", " num_cells = len( all_cells )\n", " positions = np.zeros( (num_cells,2) ); \n", " for n in range(num_cells):\n", " positions[n,:] = all_cells[n].position \n", " plt.plot( positions[:,0] , positions[:,1] , 'wo' )\n", " plt.axis('image')\n", " plt.show();\n", " plt.pause(0.0001)\n", " return; \n", "\n", "fullscreen_setup(1,0.01)\n", "\n", "# now make 10 cells with random positions \n", "number_of_cells = 10;\n", "center = [ 0.5*(E.a+E.b) , 0.5*(E.c+E.d)]\n", "for n in range(number_of_cells):\n", " c = Cell(); # create cell. \n", " r = 5 * np.random.normal();\n", " angle = 6.28318530718 * np.random.uniform(); \n", " c.position[0] = center[0] + r * np.cos(angle)\n", " c.position[1] = center[1] + r * np.sin(angle)\n", " \n", " c.equilibrium_spacing = 15 + 10*np.random.uniform(); \n", " all_cells.append(c)\n", " \n", "plot_cells( E, all_cells )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's run the mechanics steps a few times and plot." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "dt = 0.1 \n", "\n", "fullscreen_setup(1,0.001)\n", "\n", "for n in range( 250 ):\n", " for c in all_cells:\n", " c.update_velocity( dt , E, all_cells )\n", " for c in all_cells:\n", " c.update_position( dt , E, all_cells )\n", " if( n % 25 == 0 ):\n", " plot_cells( E, all_cells ) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 6. Adding cell division and death. (Cell v2) \n", "Let's add division and death to the `Cell` class. \n", "\n", "When a cell `c` divides, we want to make a fully copy of it and place it at the end of the `all_cells` data structure. \n", "\n", "Recall that in Python, the operation `b = a` doesn't actually make a new object `b` that is a copy of `a`. Instead, it makes a pointer, and changing `b` turns out to change `a`. This is stupid and deeply frustrating, but Python had good reasons to do it. We can solve it by using the `deepcopy` function in the `copy` package. \n", "\n", "When a cell `c` dies, we want to remove it from the `all_cells` data structure. \n", "\n", "Let's build these basic methods into the an upgraded `Cell` class. " ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "import copy\n", "\n", "all_cells = list(); \n", "\n", "class Cell:\n", " def __init__( self ):\n", " self.position = np.array( [0.0,0.0] ); \n", " self.velocity = np.array( [0.0,0.0] );\n", " self.mechanics_distance = 30.0; \n", " self.equilibrium_spacing = 20.0;\n", " self.mechanics_contant = 1.0; \n", " def update_velocity( self, dt , env , all_cells ):\n", " self.velocity = np.array([0.0,0.0])\n", " # spring-like interactions with cells \n", " for c in all_cells: \n", " displacement = c.position - self.position \n", " dist = np.linalg.norm( displacement )\n", " if( dist < self.mechanics_distance and dist > 1e-16 ):\n", " dv = dist - self.equilibrium_spacing \n", " displacement = displacement / dist; \n", " displacement = displacement * dv; \n", " displacement = displacement * self.mechanics_contant; \n", " self.velocity = self.velocity + displacement \n", " angle = 6.28318530718 * np.random.uniform(); \n", " perturbation_size = 0.1 * self.mechanics_contant; \n", " self.velocity[0] += perturbation_size * np.cos(angle)\n", " self.velocity[1] += perturbation_size * np.sin(angle)\n", " return\n", " def division( self, all_cells ):\n", " # make a brand new cell\n", " c = copy.deepcopy( self ); \n", " # append it to the data structure \n", " all_cells.append( c ); \n", " # set its position to be near its parent cell\n", " r = self.equilibrium_spacing; \n", " angle = np.random.uniform()\n", " c.position[0] = r*np.cos(angle)\n", " c.position[1] = r*np.sin(angle)\n", " return; \n", " def death( self, all_cells ):\n", " all_cells.remove( self )\n", " del self; \n", " def update_position( self, dt , env, all_cells ):\n", " self.position = self.position + dt*self.velocity; \n", " return; " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's test it! \n", "\n", "we'll create 10 cells. Every 50 time steps, we'll divide a random cell. Every 200 steps, we'll kill a random cell " ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "all_cells.clear() \n", "\n", "# now make 10 cells with random positions \n", "number_of_cells = 10;\n", "center = [ 0.5*(E.a+E.b) , 0.5*(E.c+E.d)]\n", "for n in range(number_of_cells):\n", " c = Cell(); # create cell. \n", " r = 5 * np.random.normal();\n", " angle = 6.28318530718 * np.random.uniform(); \n", " c.position[0] = center[0] + r * np.cos(angle)\n", " c.position[1] = center[1] + r * np.sin(angle)\n", " all_cells.append(c)\n", "\n", "fullscreen_setup(1,0.01) \n", "dt = 0.1; \n", "\n", "for n in range( 1000 ):\n", " for c in all_cells:\n", " c.update_velocity( dt , E, all_cells )\n", " for c in all_cells:\n", " c.update_position( dt , E, all_cells )\n", " if( n % 50 == 0 ):\n", " n_cells = len( all_cells )\n", " n_divide = np.random.randint( n_cells ); \n", " all_cells[n_divide].division( all_cells )\n", " if( n % 200 == 0 ):\n", " n_cells = len( all_cells )\n", " n_die = np.random.randint( n_cells ); \n", " all_cells[n_die].death( all_cells )\n", " if( n % 100 == 0 ):\n", " plot_cells( E, all_cells ) \n", " \n", " " ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "for n in range( 1000 ):\n", " for c in all_cells:\n", " c.update_velocity( dt , E, all_cells )\n", " for c in all_cells:\n", " c.update_position( dt , E, all_cells )\n", " if( n % 100 == 0 ):\n", " plot_cells( E, all_cells ) \n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 7. Update cells to probabilistically do division and death. (Cell v3) \n", "\n", "If the cell's birth rate is $b$, then the probability of the cell dividing in $[t,t+\\Delta t]$ is: \n", "\n", "$$\\textrm{Prob}_\\textrm{birth} = b \\Delta t.$$\n", "\n", "Similarly, if the death rate is $d$, then \n", "\n", "$$\\textrm{Prob}_\\textrm{death} = d \\Delta t.$$\n", "\n", "Suppose event $X$ has probability $p$. Then in numerics if you want to evaluate whether event $X$ happens, you use a uniform random number.\n", "\n", "Let $0 \\le u \\le 1$ be a uniformly distributed random number. If $u \\le p$, then $X$ happens. If $u > p$, then $X$ does not happen. \n", "\n", "Let's use this to automate birth and death in a revised `Cell` class. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "import copy\n", "\n", "all_cells = list(); \n", "\n", "class Cell:\n", " def __init__( self ):\n", " self.position = np.array( [0.0,0.0] ); \n", " self.velocity = np.array( [0.0,0.0] );\n", " self.mechanics_distance = 30.0\n", " self.equilibrium_spacing = 20.0;\n", " self.mechanics_contant = 1.0;\n", " self.birth_rate = 0.01; \n", " self.death_rate = 0.005; \n", " def update_velocity( self, dt , env , all_cells ):\n", " self.velocity = np.array([0.0,0.0])\n", " # spring-like interactions with cells \n", " for c in all_cells: \n", " displacement = c.position - self.position \n", " dist = np.linalg.norm( displacement )\n", " if( dist < self.mechanics_distance and dist > 1e-16 ):\n", " dv = dist - self.equilibrium_spacing \n", " displacement = displacement / dist; \n", " displacement = displacement * dv; \n", " displacement = displacement * self.mechanics_contant; \n", " self.velocity = self.velocity + displacement \n", " angle = 6.28318530718 * np.random.uniform(); \n", " perturbation_size = 0.1 * self.mechanics_contant; \n", " self.velocity[0] += perturbation_size * np.cos(angle)\n", " self.velocity[1] += perturbation_size * np.sin(angle)\n", " return\n", " def division( self, all_cells ):\n", " # make a brand new cell\n", " c = copy.deepcopy( self ); \n", " # append it to the data structure \n", " all_cells.append( c ); \n", " # set its position to be near its parent cell\n", " r = self.equilibrium_spacing; \n", " angle = np.random.uniform()\n", " c.position[0] = r*np.cos(angle)\n", " c.position[1] = r*np.sin(angle)\n", " return; \n", " def death( self, all_cells ):\n", " all_cells.remove( self )\n", " del self\n", " def update( self, dt, env, all_cells ):\n", " prob_birth = self.birth_rate * dt; \n", " prob_death = self.death_rate * dt; \n", " if( np.random.uniform() <= prob_birth ):\n", " self.division( all_cells );\n", " return; \n", " if( np.random.uniform() <= prob_death ):\n", " self.death( all_cells );\n", " return; \n", " return; \n", " def update_position( self, dt , env, all_cells ):\n", " self.position = self.position + dt*self.velocity; \n", " return; \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's test it! " ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "all_cells.clear() \n", "\n", "# now make 10 cells with random positions \n", "number_of_cells = 10;\n", "center = [ 0.5*(E.a+E.b) , 0.5*(E.c+E.d)]\n", "for n in range(number_of_cells):\n", " c = Cell(); # create cell. \n", " r = 5 * np.random.normal();\n", " angle = 6.28318530718 * np.random.uniform(); \n", " c.position[0] = center[0] + r * np.cos(angle)\n", " c.position[1] = center[1] + r * np.sin(angle)\n", " all_cells.append(c)\n", " \n", "fullscreen_setup(1,0.01) \n", "dt = 0.1; \n", "\n", "for n in range( 1000 ):\n", " for c in all_cells:\n", " c.update( dt , E, all_cells ); \n", " for c in all_cells:\n", " c.update_velocity( dt , E, all_cells )\n", " for c in all_cells:\n", " c.update_position( dt , E, all_cells )\n", " if( n % 50 == 0 ):\n", " plot_cells( E, all_cells ) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Section 8. Coupling cells to the environment (Cells v4) \n", "Now, we want cells to absorb substrate from the environment, and increase their birth rate when they have more substrate, and die if substrate is too low. \n", "\n", "We'll need to be able to figure out where the cell is in the environment to sample it. Here's a way to get the $(i,j)$ index of the cell\n", "\n", "$$i = \\textrm{round}\\left( \\frac{ x-a}{\\Delta x} \\right).$$\n", "$$j = \\textrm{round}\\left( \\frac{ y-c}{\\Delta y} \\right).$$\n", "\n", "Let's extend the `Cell` class to sample the substrate $c$, consume substrate from the environemnt, and change birth and death rates based on the substrates. \n", "\n", "For the birth rate, let's set $b = b_0 \\frac{c-c_\\textrm{death}}{1-c_\\textrm{death}}$ where $b_0$ is the base birth rate. \n", "\n", "And for the death rate, we'll set $d = 9e99$ if $c \\le c_\\mathrm{death}.$ \n", "\n", "For substrate consumption, let's use the (unconditionally stable) backwards Euler method: \n", "\n", "$$\\frac{ C^{n+1}_{i,j} - C^n_{i,j}}{\\Delta t} = -r_c C^{n+1}_{i,j},$$\n", " \n", "or \n", "\n", "$$C^{n+1}_{i,j} = \\frac{1}{1+r_c \\Delta t} C^n_{i,j}.$$\n", "\n", "Let's update our `Cell` class!" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "import copy\n", "\n", "all_cells = list(); \n", "\n", "class Cell:\n", " def __init__( self ):\n", " self.position = np.array( [0.0,0.0] ); \n", " self.velocity = np.array( [0.0,0.0] );\n", " self.mechanics_distance = 20.0; \n", " self.equilibrium_spacing = 15.0;\n", " self.mechanics_contant = 10.0;\n", " self.birth_rate = 0.01; \n", " self.death_rate = 0.001; \n", " self.death_threshold = 0.3; \n", " self.consumption_rate = 1.0; \n", " def update_velocity( self, dt , env , all_cells ):\n", " self.velocity = np.array([0.0,0.0])\n", " # spring-like interactions with cells \n", " for c in all_cells: \n", " displacement = c.position - self.position \n", " dist = np.linalg.norm( displacement )\n", " if( dist < self.mechanics_distance and dist > 1e-16 ):\n", " dv = dist - self.equilibrium_spacing \n", " displacement = displacement / dist; \n", " displacement = displacement * dv; \n", " displacement = displacement * self.mechanics_contant; \n", " self.velocity = self.velocity + displacement \n", " angle = 6.28318530718 * np.random.uniform(); \n", " perturbation_size = 0.1 * self.mechanics_contant; \n", " self.velocity[0] += perturbation_size * np.cos(angle)\n", " self.velocity[1] += perturbation_size * np.sin(angle)\n", " return\n", " def division( self, all_cells ):\n", " # make a brand new cell\n", " c = copy.deepcopy( self ); \n", " # append it to the data structure \n", " all_cells.append( c ); \n", " # set its position to be near its parent cell\n", " r = self.equilibrium_spacing; \n", " angle = np.random.uniform()\n", " c.position[0] = r*np.cos(angle)\n", " c.position[1] = r*np.sin(angle)\n", " return; \n", " def death( self, all_cells ):\n", " all_cells.remove( self )\n", " del self\n", " def update( self, dt, env, all_cells ):\n", " # find my coordinates in the environment\n", " i = int( np.round( (self.position[0]-env.a)/env.dx ) ); \n", " j = int( np.round( (self.position[1]-env.c)/env.dy ) ); \n", " # consume substrate (backwards Euler)\n", " constant = 1.0 + dt*self.consumption_rate; \n", " env.C[i,j] /= constant; \n", " # update birth and death rates based on C \n", " substrate = env.C[i,j]\n", " birth_rate = self.birth_rate * (substrate - self.death_threshold)/(1-self.death_threshold); \n", " death_rate = self.death_rate; \n", " if( birth_rate < 0 ):\n", " birth_rate = 0.0; \n", " if( env.C[i,j] < self.death_threshold ):\n", " self.death(all_cells)\n", " return; \n", " prob_birth = birth_rate * dt; \n", " prob_death = death_rate * dt; \n", " if( np.random.uniform() <= prob_birth ):\n", " self.division( all_cells );\n", " return; \n", " if( np.random.uniform() <= prob_death ):\n", " self.death( all_cells );\n", " return; \n", " return; \n", " def update_position( self, dt , env, all_cells ):\n", " self.position = self.position + dt*self.velocity; \n", " return; \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's make a combined plotting function. " ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "# make a plotting function \n", "def plot( env, all_cells , t ):\n", " # plt.figure(1 , figsize=(10,10) )\n", " env.plot() \n", " num_cells = len( all_cells )\n", " positions = np.zeros( (num_cells,2) ); \n", " for n in range(num_cells):\n", " positions[n,:] = all_cells[n].position \n", " plt.plot( positions[:,0],positions[:,1],'wo', markersize=15, \n", " markeredgecolor='k', alpha=0.5 )\n", " plt.axis('image')\n", " plt.title('number of cells: ' + str(num_cells) + ' time: ' + str(t) )\n", " plt.show();\n", " plt.pause(0.01)\n", " return; \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And let's put everything together into one loop! " ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length scale: 3162.2776601683795\n" ] } ], "source": [ "all_cells.clear() \n", "E.setup( 10000, 0.001 , 1 , 1 )\n", "\n", "# now make 25 cells with random positions \n", "number_of_cells = 25;\n", "center = [ 0.5*(E.a+E.b) , 0.5*(E.c+E.d)]\n", "for n in range(number_of_cells):\n", " c = Cell(); # create cell. \n", " r = 50 * np.random.normal();\n", " angle = 6.28318530718 * np.random.uniform(); \n", " c.position[0] = center[0] + r * np.cos(angle)\n", " c.position[1] = center[1] + r * np.sin(angle)\n", " c.consumption_rate = 2.0 + 0.25*np.random.normal() \n", " all_cells.append(c)\n", " \n", "dt = 0.1; \n", "fullscreen_setup()\n", "\n", "for n in range( 2000+1 ):\n", " for nn in range(10):\n", " E.update( 0.1*dt )\n", " for c in all_cells:\n", " c.update( dt , E, all_cells ); \n", " for c in all_cells:\n", " c.update_velocity( dt , E, all_cells )\n", " for c in all_cells:\n", " c.update_position( dt , E, all_cells )\n", " if( n % 200 == 0 ):\n", " plot( E, all_cells , n ) \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "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.7" } }, "nbformat": 4, "nbformat_minor": 4 }