{ "cells": [ { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import random\n", "from numpy import *\n", "import matplotlib.pyplot as plt\n", "from matplotlib.pyplot import *\n", "\n", "def celldistribution(g,R,RR,celldiameter,cellnumber): \n", " positioninR = 0\n", " validposition = 0\n", " distribution = []\n", " RRdistribution = []\n", " temp = []\n", "\n", " # Place N number of points randomly on a circle with R radius \n", " # Placing the first point\n", "\n", " ang = random.random() * 2 * np.pi\n", " rad = R * sqrt(random.random())\n", " x = rad * cos(ang)\n", " y = rad * sin(ang)\n", " distribution.append([x,y])\n", "\n", " # Placing an additional N-1 point, with the consideration that they can not be closer than m\n", "\n", " while positioninR < cellnumber-1:\n", " ang = random.random() * 2 * np.pi\n", " rad = RR * sqrt(random.random())\n", " x = rad * cos(ang)\n", " y = rad * sin(ang)\n", " for i in range (0,len(distribution),1):\n", " lenx = distribution[i][0] - x\n", " leny = distribution[i][1] - y\n", " distance = sqrt(np.power(lenx,2.) + np.power(leny,2.))\n", " if distance > celldiameter:\n", " validposition += 1\n", "\n", " if validposition == len(distribution):\n", " if rad < R:\n", " distribution.append([x,y]) # R radius cirle points\n", " positioninR += 1\n", " else:\n", " temp.append([x,y])\n", " validposition = 0\n", "\n", " RRdistribution.extend(distribution) # 2*R radius cirle points\n", " RRdistribution.extend(temp) # 2*R radius cirle points\n", " temp = []\n", "\n", " np.savetxt(\"distribution-%d.txt\" % (g),distribution)\n", " np.savetxt(\"RRdistribution-%d.txt\" % (g),RRdistribution)\n", " \n", " return distribution,RRdistribution\n", "\n", "\n", "# Calculating the distance between the cells\n", "\n", "def celldistance(g,RRdistribution): \n", " RRdist = []\n", " temp = []\n", "\n", " for i in range (0,len(RRdistribution),1): \n", " for j in range (0,len(RRdistribution),1):\n", " lenx = RRdistribution[i][0] - RRdistribution[j][0]\n", " leny = RRdistribution[i][1] - RRdistribution[j][1]\n", " distance = sqrt(np.power(lenx,2.) + np.power(leny,2.))\n", " temp.extend([distance])\n", " RRdist.extend([temp])\n", " temp = []\n", " \n", " return RRdist \n", " \n", " \n", "# Calculating the Gaussian distribution around all the cells\n", "\n", "def cellgauss(g,RRdist,RRdistribution,B):\n", " RRgau = [] \n", " temp = []\n", " \n", " for i in range (0,len(RRdistribution),1): \n", " for j in range (0,len(RRdistribution),1):\n", " gauss = (1 / (sqrt(np.pi * B))) * np.exp(-np.power(RRdist[i][j],2.) / B)\n", " temp.extend([gauss])\n", " RRgau.extend([temp])\n", " temp = []\n", "\n", " np.savetxt(\"RRgau-%d.txt\" % (g),RRgau)\n", " \n", " return RRgau\n" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [], "source": [ "def saveddata(g):\n", " RRdistribution = []\n", " distribution = []\n", " RRgau = []\n", " \n", " with open('distribution-%d.txt' % (g)) as f1:\n", " w1 = [float(x1) for x1 in next(f1).split()] # read first line\n", " distribution.append(w1)\n", " for line in f1: # read rest of lines\n", " distribution.append([float(x1) for x1 in line.split()])\n", " f1.close()\n", " \n", " with open('RRdistribution-%d.txt' % (g)) as f2:\n", " w2 = [float(x2) for x2 in next(f2).split()] # read first line\n", " RRdistribution.append(w2)\n", " for line in f2: # read rest of lines\n", " RRdistribution.append([float(x2) for x2 in line.split()])\n", " f2.close()\n", " \n", " with open('RRgau-%d.txt' % (g)) as f3:\n", " w3 = [float(x3) for x3 in next(f3).split()] # read first line\n", " RRgau.append(w3)\n", " for line in f3: # read rest of lines\n", " RRgau.append([float(x3) for x3 in line.split()])\n", " f3.close() \n", " return distribution, RRdistribution, RRgau" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "def datafile(fit):\n", " data = []\n", " with open('fit%d-2.txt' % (fit)) as f:\n", " w = [float(x) for x in next(f).split()] # read first line\n", " data.append(w)\n", " for line in f: # read rest of lines\n", " data.append([float(x) for x in line.split()])\n", " f.close()\n", " \n", "def MML(dmax,dsteps,mu,alphaR,beta,distribution,RRdistribution,RRgau,C,cellnumber):\n", " dam = []\n", " storeD = []\n", " Ci = 0\n", " \n", " for d in np.arange (0,dmax,dsteps):\n", " celldeath = 0\n", " l = d / mu\n", "\n", " LQ = np.exp(-alphaR*d-beta*np.power(d,2)) # LQ model fit\n", "\n", " for i in range (0,len(RRdistribution),1): # Give a damage number for every cell according the the Poisson distribution\n", " damage = np.random.poisson(l)\n", " dam.extend([damage])\n", " celldamage = []\n", " for i in range (0,len(dam),1): # Only those are of concern that are from the R radius circle\n", " celldamage.extend([dam[i]])\n", "\n", " for i in range (0,len(distribution),1):\n", " for j in range (0,len(RRdistribution),1):\n", " Ci += RRgau[i][j] * celldamage[j] # Calculate the percieved signal concentration for every cell\n", " if float(celldamage[i]) > Ci / C + ms: # Make the comparison between the signal and the cell's own damage\n", " celldeath += 1 # Add 1 to the cell death counter if the cell damage is higher\n", " Ci = 0\n", "\n", " survival = ((cellnumber - celldeath) / float((cellnumber))) # Calculate the survival ratio\n", "\n", " if survival < LQ: # Use the model if the SF is smaller than the one from the LQ\n", " storeD.append([d])\n", " storesurv0.append(survival)\n", " else: # Use the LQ otherwise\n", " storeD.append([d])\n", " storesurv0.append(LQ)\n", "\n", " dam = []\n", "\n", " return storeD, storesurv0\n", "\n", "def avgdev(storeD,storesurv0,gmax,fit):\n", " avg = []\n", " deviation = []\n", "\n", "\n", " for m in range (0,len(storeD),1):\n", " averageall = []\n", " summ = 0\n", " step = 0\n", " n = 0\n", " for i in range (0,len(storesurv0),1):\n", " if i == m + len(storeD) * n:\n", " averageall.append(storesurv0[i])\n", " n += 1\n", "\n", " summ = np.sum(averageall)\n", " step = summ / gmax\n", " avg.append(step)\n", "\n", " vn = 0\n", " var = 0\n", " dev = 0\n", " for i in range (0,len(averageall),1):\n", " vn += np.power(averageall[i] - step,2)\n", " var = vn / gmax\n", " dev = sqrt(var)\n", " deviation.append(dev)\n", "\n", " np.savetxt(\"Dose.txt\",storeD)\n", " np.savetxt(\"SF-average%d.txt\" % (fit),avg)\n", " np.savetxt(\"SF-deviation%d.txt\" % (fit),deviation)\n", " \n", " return avg, deviation\n", " \n", "def plot(storeD,avg,deviation):\n", "\n", " plt.clf()\n", " plt.cla()\n", " plt.plot(storeD,avg,'b.') # Model data plot\n", "# z, w = np.loadtxt('fit%d.txt' % (fit), unpack=True) # Fit data\n", "# plt.plot(z,w,'ro',label='Data') # Fit data plot\n", " plt.errorbar(storeD, avg, deviation, linestyle='None', marker='^')\n", "\n", " xlabel (\"D [Gy]\"); ylabel (\"Surviving fraction\") # Labels\n", " plt.title (\"Surviving fraction\") # Title\n", " plt.grid(True) # Grid\n", " plt.show()\n", " plt.clf()\n", " plt.cla()\n", " plt.close()" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "\"\"\"Generate the system setup\"\"\"\n", "\n", "\"\"\"Fixed parameters\"\"\"\n", "\n", "R = 0.5 # Size of the system\n", "RR = 1 # 2*R size\n", "B = 0.02 # Parameter defining the half width of the Poisson distribution\n", "celldiameter = 0.006 # Size of the cell nucleus\n", "cellnumber = 600 # Number of cells\n", "C = 196 # Normalization factor for the cell signal\n", "\n", "\"\"\"Number of simulations\"\"\"\n", "\n", "g = 0\n", "gmax = 10 # Calculations g times\n", "\n", "datafile(fit)\n", "\n", "while g < gmax:\n", " distribution, RRdistribution = celldistribution(g,R,RR,celldiameter,cellnumber)\n", " RRdist = celldistance(g,RRdistribution)\n", " RRgau = cellgauss(g,RRdist,RRdistribution,B)\n", " g += 1" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwcZb3v8c+XBBFIBBGJaFAQcAFUJBHEjQl6laiI+vLI4gHlyoEcxe0cuLijR1Be6vEosoToReQKDIoiqFHgQEaQTUkOW2RJ2EKMGjAMSWBCmMnv/lE1nU7T013T09XVy/f9es2L6eqq6t80lf718/yeeh5FBGZmZgCbFR2AmZm1DycFMzMrcVIwM7MSJwUzMytxUjAzsxInBTMzK3FSsJ4gaa6kL2XY78WS1kqalEMMW0r6laTHJf2s2eev8nqZ/mazcvJ9ClYUSW8CvgnsCYwAdwGfjog/FRpYTiQdCXwCeENEDDf53B8BjomINzXzvNZ7JhcdgPUmSc8Bfg38K/BT4FnAm4GnGjiXSL7gbGhqkM33EuDesRKCpMnNThZm4+XuIyvKywAi4qKIGImIoYi4MiJuB5D0FUk/Gd1Z0s6SQtLk9PGApFMlXQ88CXxe0i3lLyDpM5IuT38/T9Ip6e93SXp32X6TJT0qaZ8xXudrkq6XtEbSlZK2Lzv2KEkPSfqHpC9JelDS2yr/WElfBb4MHJp2T31U0kfS8/6XpFXAVyTtKuma9HyPSrpA0rZl59lJ0i8kPZLuc4akVwJzgf3Tcw9W/s3p43+RtFTSKkmXS3ph2XMhaY6kJZIek3RmmmytxzgpWFHuBUYk/VjSbEnPbeAcRwLHAlOB7wMvl7R72fNHABdWOe4i4PCyx+8AHo2IRWO8zhHA0cAOJC2aEwAk7QGcBXwI2BHYBnhRtRNExMnA14GLI2JKRPzf9Kn9gPvTc58KCPgG8ELglcBOwFfS15tE0rp6CNg5fa3+iLgLmAPcmJ67lERGSTowPe8H01gfAvordns38DrgNel+7xjj/bAu5qRghYiI1cCbgAB+ADySfnudNo7TnBcRiyNiOCIeBy4j/bBPk8MrgMurHHch8B5JW6WPx0oeo34UEfdGxBBJV9fe6fYPAL+KiD9ExHqSlsB4i3QrIuL76d8wFBFLI+KqiHgqIh4BvgMckO67L0myODEinoiIdRHxh4yv8yHg3IhYFBFPAZ8jaVnsXLbPaRExGBHLgAVlf6f1ECcFK0xE3BURH4mI6cBeJB943x3HKR6ueHwhG1sARwC/jIgnq7zuUpKi9sFpYngPtZPC38p+fxKYkv7+wvIY0tf6xzjih4q/QdIOkvol/UXSauAnwGh31U7AQw3WHV5I0joYjXVtGmt5y2asv9N6iJOCtYWIuBs4jyQ5ADwBbFW2ywuqHVbx+Epge0l7kySHWh/0o11IhwB/ThPFeP0VmD76QNKWwPPGeY7Kv+Eb6bZXR8RzgH8m6VKCJIG8eLTeUec8lVaQFLpHY906jfUv44zXupyTghVC0isk/buk6enjnUg+pG9Kd7kVeEt638A2JN0dNaXfoC8BvgVsB1xVY/d+4O0ko59qJY9aLiFpbbxB0rOAr7LxA7xRU4G1wKCkFwEnlj33R5JEdJqkrSU9W9Ib0+f+DkxP46jmQuBoSXtL2oKkvnFzRDw4wXityzgpWFHWkBRZb5b0BEkyuBP4d4CIuAq4GLgdWEhSYM3iQuBtwM9qdbNExF+BG4E3pK8zbhGxmOS+g36SD+s1wEoaGFZb5qvAPsDjwG+AX5S93ghwMLAbsAxYDhyaPn0NsBj4m6RHq8R6NfAl4OdprLsCh00gTutSvnnNrEkkTQEGgd0j4oGi4zFrhFsKZhMg6WBJW6V99N8G7gAeLDYqs8Y5KZhNzCEkRdwVwO7AYeHmt3Uwdx+ZmVmJWwpmZlbScRPibb/99rHzzjs3dOwTTzzB1ltv3dyAmsSxNaadY4P2js+xNaZTY1u4cOGjEfH8uieJiI76mTFjRjRqwYIFDR+bN8fWmHaOLaK943NsjenU2IBbIsNnrLuPzMysxEnBzMxKnBTMzKzEScHMzEqcFMzMrCS3pCDpXEkrJd05xvOSdHq6PODtkvbJKxaAlavX8fWbh1i5Zl2eL2Nm1tHybCmcBxxU4/nZJNMC7E6ypOLZOcbC5y9cwr2rNvCFCxqZNt/MrDfklhQi4lpgVY1dDgHOT4fQ3gRsK2nHPGKZv2AdVy1dDoIrlzzMbxe4tWBmVk2ucx+l67/+OiL2qvLcr0nWhP1D+vhq4KSIuKXKvseStCaYNm3ajP7+yvXGa/viLzfn4cmPo8lBDIudhrfhlPc+Pe6/J09r165lypT2XP3QsTWuneNzbI3p1NhmzZq1MCJm1jtHkdNcVFuhqmqGioh5wDyAmTNnRl9fX+YXWbl6HX+/cgHakJxak4OVz1rNHjNmscPUZ4876LwMDAwwnr+rlRxb49o5PsfWmG6PrcjRR8tJFiIfNZ1k+uGmOv3qJaBNc00oOP1q1xbMzCoVmRQuB45KRyG9Hng8kiUSm2rRskGeHtk0KTw9Eix66LFmv5SZWcfLrftI0kVAH7C9pOXAycDmABExF5gPvBNYCjwJHJ1HHPM/9ebS7+3c7DMzawe5JYWIOLzO8wF8PK/XNzOz8fMdzWZmVuKkYGZmJU4KZmZW4qRgZmYlTgpmZlbipGBmZiVOCmZmVuKkYGZmJU4KZmZW4qRgZmYlTgpmZlbipGBmZiVOCmZmVuKkkHrfmdfzqpOvYOUar99sZr3LSSF137Jh1qwb5gsXeEU2M+tdTgrA/AXreHzkSRBcueRhfrvArQUz601OCsAZC5YA6ZKdCr6/wK0FM+tNPZ8UVq5ex5Lh5WhykhQ0OVg6/LBrC2bWk3o+KZx+9RJQbLItFJx+tVsLZtZ7ej4pLFo2yNMjmyaFp0eCRQ89VlBEZmbFmVx0AEWb/6k3Fx2CmVnb6PmWgpmZbeSkYGZmJU4KZmZW4qRgZmYlTgpmZlbipGBmZiVOCmZmVuKkYGZmJU4KZmZW4qRgZmYlTgpmZlbipGBmZiVOCmZmVuKkYGZmJU4KGb3vzOt51clXeEU2M+tquSYFSQdJukfSUkmfrfL8cyVdKul2SX+UtFee8UzEfcuGWbNumC9c4BXZzKx75ZYUJE0CzgRmA3sAh0vao2K3zwO3RsSrgaOA7+UVz0TMX7COx0eeBMGVSx7mtwvcWjCz7pRnS2FfYGlE3B8R64F+4JCKffYArgaIiLuBnSVNyzGmhpyxYAmQLtmp4PsL3Fows+6kiKi/VyMnlj4AHBQRx6SPjwT2i4jjy/b5OvDsiPg3SfsCN6T7LKw417HAsQDTpk2b0d/f31BMa9euZcqUKeM6ZnDdBk4YGGK4bNtkwbf7tmTbLZqXUxuJrVUcW+PaOT7H1phOjW3WrFkLI2JmvXPkuUazqmyrzECnAd+TdCtwB/A/sMnnb3JQxDxgHsDMmTOjr6+voYAGBgYY77FfvPQONOlhGNkYujYTtwztwCnvaF4JpJHYWsWxNa6d43Nsjen22PJMCsuBncoeTwdWlO8QEauBowEkCXgg/Wkbi5YN8vTIprns6ZFg0UOPFRSRmVl+8kwKfwJ2l7QL8BfgMOCI8h0kbQs8mdYcjgGuTRNF25j/qTcXHYKZWcvklhQiYljS8cAVwCTg3IhYLGlO+vxc4JXA+ZJGgD8DH80rHjMzqy/PlgIRMR+YX7FtbtnvNwK75xmDmZll5zuazcysxEnBzMxKnBTMzKzEScHMzEqcFMzMrMRJwczMSuoOSZX0MuBE4CXl+0fEgTnG1XEOPedGAC4+bv+CIzEza1yW+xR+BswFfgCM5BuOmZkVKUtSGI6Is3OPxMzMCpelpvArSR+TtKOk7UZ/co+sw6wf3sCfV6z2cp1m1tGyJIUPk9QUbgAWpj+35BlUJ/JynWbWDep2H0XELq0IpJONLtepyaPLde7G7FnPLjosM7Nxq9tSkLS5pE9KuiT9OV7S5q0IrlN4uU4z6xZZuo/OBmYAZ6U/M9JtBqxcvY4lw8vR5CQpaHKwdPhh1xZSh55zY2m4biPPm1lrZUkKr4uID0fENenP0cDr8g6sU5x+9RLQpiuzhYLTr3ZrIQsX6M3aS5akMCJp19EHkl6K71co8XKdE+MCvVl7yXKfwonAAkn3AyK5s/noXKPqIKPLdfqO5vFzgd6s/WQZfXS1pN2Bl5Mkhbsj4qncI+swTgbjV61AP3vWXoXGZNbrxkwKkg6MiGskvb/iqV0lERG/yDk262JjF+h3Y4epbi2YFaVWTeGA9L8HV/l5d85xWZcYq5CctUDv0UlmrTVmSyEiTk5//Y+IeKD8OUm+oc2A+rWU+5YNs4akkPyDORu7hlygN2tPWQrNPwf2qdh2Ccn9CmZjqlVIzlqgXz+8gaUr17JyzTp3K5m1QK2awiuAPYFtKuoKzwH8r9PqylJIrlegH6ulYWb5qFVTeDlJ7WBbNq0n7AP8S/6hWSdrxp3eoy0NNNrS8A1uZnmrVVO4DLhM0v4R4UqfjUutQvIp7832jb9eS+PQc25kcHCIvr4mBW1mme5oniNp29EHkp4r6dwcY7IuMNFCcpaWxvrhDSxbvcFTZJg1UZZC86sjYnD0QUQ8Jum1OcbUlXrtjueJ3umdpaVx37JhhsD1BrMmytJS2EzSc0cfpKuuZUkm1gPymtCuXkvD9QazfGT5cP9P4AZJl6SP/wk4Nb+QrJPkNTqoXkvDU2SY5SPL3EfnS1oIzCKZ++j9EfHn3COztteKCe2qdTt5igyz/GTpPiIiFgM/BS4D1kp6ca5RdaFuXDcgy4pzFx+3f9PrKF7Dwiw/WZbjfI+kJcADwO+BB4Hf5hxX1+m2dQOKXHHOU2SY5SdLTeFrwOuB/46I10qaBRyeb1jdpRvXDWjGfQiNKq83DA4OcsVJs6vu12sjvsyaIUv30dMR8Q+SUUibRcQCYO+c4+oqWbpZOk07fFu/+Lj9+dx+W7bs9cx6QZaWwqCkKcC1wAWSVgLD+YbVPbq1KOoV58y6U5aWwiHAk8BngN8B95HMgWQZuChanG4s7pvlrWZSkDQJuCwiNkTEcET8OCJOT7uT6pJ0kKR7JC2V9Nkqz28j6VeSbpO0WFLXrf3cDt0svarbivtmrVCz+ygiRiQ9KWmbiHh8PCdOE8qZwP8ClgN/knR5xT0OHwf+HBEHS3o+cI+kCyJi/Tj/jrbV6d0snTrpXDcW981aIUtNYR1wh6SrgCdGN0bEJ+scty+wNCLuB5DUT9IVVZ4UApgqScAUYBWuV1gT+I5ns8YoImrvIH242vaI+HGd4z4AHBQRx6SPjwT2i4jjy/aZClwOvAKYChwaEb+pcq5jgWMBpk2bNqO/v79mzGNZu3YtU6ZMaejYvLVrbN+4eYiRkRG++IbqsX3j5iGAwkYBVXvfBtdt4ISBoU2+XUwWfLtvS7bdYmOPaStib9f/r+DYGtWpsc2aNWthRMysd45aK69dHRFvBfaIiJMaiE9VtlVmoHcAtwIHArsCV0m6LiJWb3JQxDxgHsDMmTOjr8G+jIGBARo9Nm/tGtvZ9yT3AowVW9EhV3vfvnjpHWjSw1BWy9Fm4pahHTjlHRtbC2ffk3Tp9fXl16XXrv9fwbE1qttjq9V9tKOkA4D3pF0/m3zIR8SiOudeDuxU9ng6sKJin6OB0yJpriyV9ABJq+GPWYI3q8bFfbPG1UoKXwY+S/Jh/p2K54Lk230tfwJ2l7QL8BfgMOCIin2WAW8FrpM0jWQJ0Puzhd49OrWY2646vbhvVqRay3FeAlwi6UsR8bXxnjgihiUdD1wBTALOjYjFkuakz88lmULjPEl3kLREToqIRxv5Q8yazUnFelGWqbPHnRDKjp0PzK/YNrfs9xXA2xs9v1kt9T7M1w9vYOnKtaxcs66j7y43a6ZMU2ebdSPf3Gb2TE4K1pO8nKdZdVnWU9iuys/mrQjOLC9ZZq713EnWi7K0FBYBjwD3AkvS3x+QtEjSjDyDM8tD1gWC3L1kvShLUvgd8M6I2D4ingfMJlma82PAWXkG1yvWD29g2eoN/kbaIllmrnX3kvWqLElhZkRcMfogIq4E3hIRNwFb5BZZD7lv2TBDw/gbaYtkubmtGxdGMssiy4R4qySdBIxOOHQo8Fg6C+qG3CLrEe0+m2d5K6Zbhm3Wu7mtWxdGMssiS0vhCJK7mn8JXAa8ON02CfhgfqH1hnb/RtqLrRgvjGS9LMvNa48Cnxjjaf8rmYB2/0ba7q2YvHjuJOtldZOCpJcBJwA7l+8fEfXmPrI6an0jPeW9xc/936trEnjuJOtlWWoKPwPmAj8ERvINp7e08zfSdm/FmFk+siSF4Yg4O/dIelD5N9LBwUGuOGl21f3y/sZa7fzt3oppBrcAzJ4pS1L4laSPAZcCT41ujIhVuUVlhWvnVkyr1Esa7l6ybpQlKYwux3li2bYAXtr8cKxdZG3F9DLPsmrdqO6Q1IjYpcqPE4L1PE+DYd2o1hrNB0bENZLeX+35iPhFfmGZtbcsw3W9op51olrdRwcA1wAHV3kuACeFFnE3Rfvp1eG61v1qLcd5cvrrMRHhoagFum/ZMGtIuil+MMcfPEXzcF3rZlmmuXhA0jxJb5Wk3CPqQRcftz+f22/Lqs95ts7242kwrJtlSQovB/4b+DhJgjhD0pvyDctGtfvcSL3Iw3Wtm2WZ+2iIZP2En0p6LvA94PckE+JZjlrVTeGaxfiMDtetpxtnmLXul2mNZkkHSDqLZBW2Z+PZUVuiVd0UHlqZj16cYdY6X5Y1mh8APg1cB+wVER+MiJ/nHpm1pJvCNYt81HtfDz3nxtId0WbtJMsdza+JiNW5R2LP0IrZOj20Mh9+X61T1bp57f9ExDeBU6WKPgwgIj6Za2SWOw+tzIffV+tktbqP7kr/ewuwsMqPdTgPrcyH31frZLVuXvtV+uvtEfE/LYrHWihLzeLi4/ZnYGCgxZF1tizvq0d8WbvKUlP4jqQdSRbb6Y+IxTnHZBXymprZK4zlI8sMs75L3dpVlllSZwF9wCPAPEl3SPpi3oFZ83ikS3vxiC9rZ5nuU4iIv0XE6cAc4Fbgy7lGZePiD/3O4rvUrZ1luU/hlZK+IulO4AzgBmB67pGZdaGxRya5tWDtIUtL4UfAY8DbI+KAiDg7IlbmHJdZx6s20aFHJlm7q1loljQJuC8ivteieMy6mifTs3ZXMylExIik50l6VkSsb1VQNj4e3tg5so748ogwK0qWIakPAddLuhx4YnRjRHwnt6hsXOoNb6yXNPzBY2ajstQUVgC/TvedWvZjbSDL8EbPgtp51g9v4M8rVrsAbS2XZT2FrzZ6ckkHkay/MAn4YUScVvH8icCHymJ5JfD8iFjV6Gv2mnoTr2VZYN5ar17rzDe3WVGyDEldIOmayp8Mx00CzgRmA3sAh0vao3yfiPhWROwdEXsDnwN+74SQXZbhjR4T33l8c5sVKUv30QnAienPl0huXrslw3H7Aksj4v60SN0PHFJj/8OBizKc11L1hjd6THxnypLIfcOi5UURz5gVu/5B0u8j4oA6+3wAOCgijkkfHwnsFxHHV9l3K2A5sFu1loKkY4FjAaZNmzajv79/3DEDrF27lilTpjR0bN4aie3L1z/JsjXP/P/34qniP964FT9e/BTXLh+mfATkJMEB0ydz1J5b5Bpbq7RzbDD++AbXbeCEgSGGy7ZNFny7b0u23WLjd7hv3DwE8Iz7IPKMrZUcW2NqxTZr1qyFETGz3jnq1hQkbVf2cDNgBvCCDPGpyraxMtDBwPVjdR1FxDxgHsDMmTOjr68vw8s/08DAAI0em7dGYrs23X2s4YvfvO06RirWRxoJ+NvwVvT1ZVtnuNHYWqWdY4Pxx/fFS+9Akx6mPJNrM3HL0A6c8o6NtYWz70n+n/f1NT5yrJ3fO8fWmGbElmVI6kKSD3MBw8ADwEczHLcc2Kns8XSSkUzVHIa7jprOs6B2Ht/cZkXLMvpolwbP/Sdgd0m7AH8h+eA/onInSdsABwD/3ODrmHWN0URej29YtLyMWWiW9DpJLyh7fJSkyySdXtGlVFVEDAPHA1eQrOL204hYLGmOpDllu74PuDIinqh2Hqvv4uP2dyugx/jeE8tLrdFH5wDrASS9BTgNOB94nLR/v56ImB8RL4uIXSPi1HTb3IiYW7bPeRFxWKN/gFmvyTJk1aOTrFG1ksKkssLvocC8iPh5RHwJ2C3/0MysGt97YnmqVVOYJGly2g30VtIhoRmOszbjrqXuMfa9J7u5tmBNUaulcBHwe0mXAUPAdQCSdiPpQjKzFvN6DJa3Mb/xR8Spkq4GdiQpBI9eiZsBn2hFcGa2KQ9ZtbzVW0/hpirb7s0vHDOrJeu9J/WGrB56zo0MDg5R7T6nid7X4vtiOluWuY/MrMO085DV9515Pa86+QrPwdWmnBTMukzRs6zWWwuinROWOSmYdaRaNywWPWS11od+0QnL6nNSMOsizZguvd43/Vo3xtX70M+SsFauXsfXbx5y91JBnBTMukjWIavrhzewbPWGqh+8E+neqfWhnzVhff7CJdy7aoO7lwripGDWRbIOWb1v2TBDwzzjg3ci3Tv1PvSzJKz5C9Zx1dLl7l4qkJOCWReZ/6k38+Bp72K/XbZjv12248HT3sWDp71rk9lXa33wT6QeUe9DP0vCKroeYp6uwqznVPvgnT1rrwlPoVHvQ7/etOCewqM9uKVg1kNqdfFMdAqNLK2UWrK8/srV6/jgOTe6CJ0jJwWzHlLrgzdrPaLe6KRGZXn9z1+4hD8+sMpF6By5+8isC411D0OtD96sU2jct2yYNSSjk34wZ69nPN/o9BblLYpqaw2PFqE1ebQWshuzZ7lbqdmcFMx6SPkH/+DgIFecNHt8x6dF6iI+mMeqhVhzufvIzDIranRQM27Ks2ycFMwskyI/mLMWwV2InjgnBTPbxFjzKhW5wE/WIrgL0RPnmoKZZVLkAj9ZhrW6EN0cTgpmPeji4/ZnYGBgXMdkHZ1UFBeim8PdR2bW8VyIbh4nBTPreEXWO7qNk4KZdbys9Q6PTqrPNQUzG5d2qyVAtkI0pKOTHlw15t3Y5paCmfWI0loNeK2GWpwUzKwneK2GbJwUzKzreXRSdk4KZtb1PE1Gdk4KZtb1PE1Gdh59ZGZdz9NkZOeWgpkZLkSPclIws57nQvRGTgpm1vPGU4j++s1DXZ0snBTMrOeNpxB976oNXV2IzrXQLOkg4HvAJOCHEXFalX36gO8CmwOPRsQBecZkZlbJheiNcmspSJoEnAnMBvYADpe0R8U+2wJnAe+JiD2Bf8orHjOzieiVQnSe3Uf7Aksj4v6IWA/0A4dU7HME8IuIWAYQEStzjMfMrCG9VIhWRNTfq5ETSx8ADoqIY9LHRwL7RcTxZfuMdhvtCUwFvhcR51c517HAsQDTpk2b0d/f31BMa9euZcqUKQ0dmzfH1ph2jg3aOz7Hlt2PFz/FtcuHKS87TBIcMH0yR+25BQCD6zZw1m1P8bG9t2DbLYop19Z632bNmrUwImbWO0eeNQVV2VaZgSYDM4C3AlsCN0q6KSLu3eSgiHnAPICZM2dGX19fQwENDAzQ6LF5c2yNaefYoL3jc2zZffO26xiJ1ZtsGwn42/BW9PUl9Yhj5t7BvY8t4zcP7FDYtNzNeN/yTArLgZ3KHk8HVlTZ59GIeAJ4QtK1wGuAezEzaxPlhehqH7zdVITOs43zJ2B3SbtIehZwGHB5xT6XAW+WNFnSVsB+wF05xmRm1nTdVITOraUQEcOSjgeuIBmSem5ELJY0J31+bkTcJel3wO3ABpJhq3fmFZOZWbONXYTejR2mdl5rIddqSETMj4iXRcSuEXFqum1uRMwt2+dbEbFHROwVEd/NMx4zs2brtmm5fUezmdkEdNu03J4628xsArrtbmi3FMzMctZJhWgnBTOzHHXa3dBOCmZmOeq0QrSTgplZjjqtEO1Cs5lZjjqtEO2WgplZwdqpEO2kYGZWoHYrRDspmJkVqN0K0U4KZmYFardCtAvNZmYFardCtFsKZmZtrpWFaCcFM7M21upCtJOCmVkby1qIbhYnBTOzNpa1EN0sLjSbmbWxLIXoZnJLwczMSpwUzMysxEnBzMxKnBTMzKzEScHMzEoUEfX3aiOSHgEeavDw7YFHmxhOMzm2xrRzbNDe8Tm2xnRqbC+JiOfXO0HHJYWJkHRLRMwsOo5qHFtj2jk2aO/4HFtjuj02dx+ZmVmJk4KZmZX0WlKYV3QANTi2xrRzbNDe8Tm2xnR1bD1VUzAzs9p6raVgZmY1OCmYmVlJ1yQFSQdJukfSUkmfrfK8JJ2ePn+7pH2yHtuC2D6UxnS7pBskvabsuQcl3SHpVkm3FBBbn6TH09e/VdKXsx7bgthOLIvrTkkjkrZLn8vtfZN0rqSVku4c4/kir7V6sRV5rdWLrbBrLWN8RV1vO0laIOkuSYslfarKPs275iKi43+AScB9wEuBZwG3AXtU7PNO4LeAgNcDN2c9tgWxvQF4bvr77NHY0scPAtsX+L71Ab9u5Ni8Y6vY/2Dgmha9b28B9gHuHOP5Qq61jLEVcq1ljK2Qay1rfAVebzsC+6S/TwXuzfPzrVtaCvsCSyPi/ohYD/QDh1TscwhwfiRuAraVtGPGY3ONLSJuiIjRFTNuAqY38fUnFFtOx+Zx/sOBi5r4+mOKiGuBVTV2KepaqxtbgddalvdtLLm/bzDu+Fp5vf01Ihalv68B7gJeVLFb0665bkkKLwIeLnu8nGe+aWPtk+XYvGMr91GSjD8qgCslLZR0bBPjGk9s+0u6TdJvJe05zmPzjg1JWwEHAT8v25zn+1ZPUdfaeLXyWsuqiGttXIq83iTtDLwWuLniqaZdc92y8pqqbKscazvWPlmOnYjM55c0i+Qf6pvKNr8xIlZI2gG4StLd6TeaVsW2iGTOlLWS3gn8Etg947F5xzbqYOD6iCj/lpfn+1ZPUddaZgVca1kUda2NVyHXm6QpJIno0xGxuvLpKoc0dM11S0thObBT2eXW8H8AAAM+SURBVOPpwIqM+2Q5Nu/YkPRq4IfAIRHxj9HtEbEi/e9K4FKS5mDLYouI1RGxNv19PrC5pO2zHJt3bGUOo6Ipn/P7Vk9R11omBV1rdRV4rY1Xy683SZuTJIQLIuIXVXZp3jWXR2Gk1T8kLZ77gV3YWEzZs2Kfd7FpIeaPWY9tQWwvBpYCb6jYvjUwtez3G4CDWhzbC9h4k+O+wLL0PSz8fUv324akH3jrVr1v6Xl3ZuyCaSHXWsbYCrnWMsZWyLWWNb6irrf0PTgf+G6NfZp2zXVF91FEDEs6HriCpNp+bkQsljQnfX4uMJ+kQr8UeBI4utaxLY7ty8DzgLMkAQxHMtPhNODSdNtk4MKI+F2LY/sA8K+ShoEh4LBIrrZ2eN8A3gdcGRFPlB2e6/sm6SKSkTLbS1oOnAxsXhZXIddaxtgKudYyxlbItTaO+KCA6w14I3AkcIekW9NtnydJ8E2/5jzNhZmZlXRLTcHMzJrAScHMzEqcFMzMrMRJwczMSpwUzMysxEnBrEw68+Wt6WyUt0n6N0nP+HdSNqPn/LJtu0v6taT70ukOFkh6S53X+4ykZZLOyOPvMRuvrrhPwayJhiJib4B0yoILSW5YOrnKvtdFxLvTfZ8N/AY4ISIuT7ftBcwExpzuICL+S9Jj6X5mhXNLwWwMkUxZcCxwvNI7k2r4EHDjaEJIj78zIs6TtJmkJZKeD5A+XppO4WDWVpwUzGqIiPtJ/p3sUGfXPUkmdKt2jg3AT0gSB8DbgNsi4tFmxWnWLE4KZvXVayU88wDp0nR1rtHJy84Fjkp//9/Aj5oVnFkzOSmY1SDppcAIsLLOrotJVu0CICLeB3wE2C59/DDwd0kHAvux6ToGZm3DScFsDGkNYC5wRtSfJOxC4I2S3lO2bauKfX5I0o3004gYaV6kZs3j0Udmm9oynYlyc2AY+H/Ad+odFBFDkt4NfEfSd4G/A2uAU8p2u5yk28hdR9a2nBTMykTEpAkcezfJ9MVjeQ1JgfnuRl/DLG/uPjJrzHpgr/Kb12qR9FmSlbM+V7H9M+m2yuUVzQrh9RTMzKzELQUzMytxUjAzsxInBTMzK3FSMDOzEicFMzMr+f8lWthspJl9HQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "\"\"\"Run the MML model\"\"\"\n", "\n", "\"\"\"Data\"\"\"\n", "\n", "fit = 68 # Data identification number\n", "ms = 0.5 # Spontanious mutation rate\n", "mu = 1 # Mutation induction rate\n", "\n", "\"\"\"LQ model parameters\"\"\"\n", "\n", "alphaR = 0.2\n", "beta = 0.05\n", "\n", "\"\"\"Dose and steps\"\"\"\n", "\n", "dmax = 2 # Maximum simulated dose\n", "dsteps = 0.05 # Simulated dose steps\n", "\n", "\"\"\"Number of simulations\"\"\"\n", "\n", "g = 0\n", "gmax = 10 # Calculations g times\n", "generated = 1 # 1 if the distribution has already been generated, 0 if not\n", "\n", "storesurv0 = []\n", "while g < gmax:\n", " if generated == 1:\n", " distribution, RRdistribution, RRgau = saveddata(g)\n", " storeD, storesurv0 = MML(dmax,dsteps,mu,alphaR,beta,distribution,RRdistribution,RRgau,C,cellnumber)\n", " g += 1\n", "avg, deviation = avgdev(storeD,storesurv0,gmax,fit)\n", "plot(storeD,avg,deviation)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }