{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Competing Yeast Species\n",
    "\n",
    "#A worksheet to facilitate analyzing the competing yeast species project in Section 7.5.3 of the text. The goal is to use the data of Table 7.2 to estimate the\n",
    "#parameters r1, K1, r2, K2, a, and b in equations (7.3)-(7.4).\n",
    "\n",
    "#Load in the data. Missing population data is indicated with a -1.\n",
    "load('comp_yeast_data.sage');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's form amalgamated data sets for runs 1 and 2, for the Saccharomyces population alone (column 2), the \n",
    "#Schizosaccharomyces population alone (column 4), the Saccharomyces population in competition (column 3), and the \n",
    "#Schizosaccharomyces population in competition (column 5).\n",
    "\n",
    "N0 = N1 + N2; #Total number of times/data points.\n",
    "times = times1+times2; #Concatenate all times\n",
    "sacc = sacc1+sacc2; #Concatenate Sacch data\n",
    "saccind = list(zip(times,sacc)); #Arrange in (time, population) pairs.\n",
    "schiz = schiz1+schiz2; #Concatenate Schiz data\n",
    "schizind = list(zip(times,schiz)); #Arrange in (time, population) pairs.\n",
    "saccm = saccmix1+saccmix2; #Concatenate Sacch data\n",
    "saccmix = list(zip(times,saccm)); #Arrange in (time, population) pairs.\n",
    "schizm = schizmix1+schizmix2; #Concatenate Schiz data\n",
    "schizmix = list(zip(times,schizm)); #Arrange in (time, population) pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We can plot the data. For example, for the saccharomyces population growing in isolation. Also, limit view to 0 to 60 in \n",
    "#time, 0 to 15 in population.\n",
    "list_plot(saccind).show(xmin=0,xmax=60,ymin=0,ymax=15);\n",
    "\n",
    "#A similar plot can be constructed for the schizosaccharomyces population in isolation, or either species in competition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Estimate Growth Parameters for Saccharomyces: We can estimate the parameters r1 and K1 for the logistic growth model for \n",
    "#the saccharomyces population using the corresponding data.\n",
    "\n",
    "#The solution to the logistic ODE u' = r1*u*(1-u/K1) with initial data u(t0) = u0 is\n",
    "var('t r1 K1 u0 t0'); #Declare t as symbolic variable, t0,u0,r1,K1 too.\n",
    "usol(t) = K1*u0/(exp(-r1*(t-t0))*(K1-u0)+u0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The initial data could be taken as any (time, population) data point, but\n",
    "#we'll use the first:\n",
    "t0 = times[0]\n",
    "u0 = sacc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The solution with this initial data is\n",
    "u1sol(t) =  K1*u0/(exp(-r1*(t-t0))*(K1-u0)+u0);\n",
    "u1sol(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot or even a glance at the data suggests that K1 is around 12 to 14. To estimate K1 and r1 more carefully,  we can use the \"find_fit\" command, that performs\n",
    "#a least squares fit of u1sol(t) to the data by adjusting r1 and K1.\n",
    "\n",
    "#But first, we must filter out the \"no data\" points by excluding those with population = -1\n",
    "saccpos =  [[times[i],sacc[i]] for i in range(N0) if sacc[i]>0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now fit the model\n",
    "sol = find_fit(saccpos,u1sol,parameters=[r1,K1],initial_guess=[0.5,14.0]) #Fit the model by adjusting r1 and K1\n",
    "r1best = sol[0].rhs();\n",
    "K1best = sol[1].rhs();\n",
    "bestu1(t) = u1sol(r1=r1best,K1=K1best) #Define the best fit function of the form in \"model\"\n",
    "print(\"Optimal r1 = \",r1best,\" Optimal K1 = \",K1best)\n",
    "print(\"Best fit u1(t) = \", bestu1(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the best fit u1(t) and the data.\n",
    "p1 = list_plot(saccind);\n",
    "p2 = plot(bestu1(t),(t,0,60),color='red');\n",
    "(p1+p2).show(xmin=0,xmax=60,ymin=0,ymax=15);\n",
    "\n",
    "#A similar procedure can be used to estimate r2 and K2 for the schizosaccharomyces data, using (for example) the data \n",
    "#point t0 = 16, u0 = 1 as an \"initial\" data point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Competition Parameters a and b: The ODEs that govern the competing species system are\n",
    "var('a b r2 K2');\n",
    "u1 = function('u1')(t);\n",
    "u2 = function('u2')(t);\n",
    "de1 = diff(u1,t) == r1*u1*(1-u1/K1 - a*u2/K1)\n",
    "de2 = diff(u2,t) == r2*u2*(1-u2/K2 - b*u1/K2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A closed-form solution cannot be obtained, but a guess-and-plot approach based on solving these ODEs numerically can work, \n",
    "#using the values for r1, K1, r2, K2 found above."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}