Predict2.0

Manifest.toml

@@ -114,6 +114,12 @@ git-tree-sha1 = "6b19601c0e98de3a8964ed33ad73e130c7165b1d"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 version = "0.22.4"
 
+[[ExponentialUtilities]]
+deps = ["LinearAlgebra", "Printf", "SparseArrays"]
+git-tree-sha1 = "1672dedeacaab85345fd359ad56dde8fb5d48a45"
+uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
+version = "1.6.0"
+
 [[FFMPEG]]
 deps = ["BinaryProvider", "Libdl"]
 git-tree-sha1 = "9143266ba77d3313a4cf61d8333a1970e8c5d8b6"
@@ -241,7 +247,7 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 version = "0.3.11"
 
 [[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Test", "UUIDs"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[PlotThemes]]

Project.toml

@@ -8,6 +8,7 @@ BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

doxorubicin.ipynb

@@ -15,7 +15,7 @@
     "conc_l, pop, g2, g1 = setup_data(\"doxorubicin\")\n",
     "\n",
     "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 4E4)\n",
-    "ps = [134.546, 0.100021, 0.0385798, 0.814924, 2.93535, 0.229024, 0.0875672, 0.168648, 0.631369, 2.37796, 88.8261, 4.03489, 2.42543]\n",
+    "ps = [144.93, 0.520304, 0.170531, 0.818625, 0.328673, 5.83542e-6, 0.123382, 0.208883, 0.638467, 10.393, 10.0785, 6.38181, 3.40993]\n",
     "effects = getODEparams(ps, conc_l)"
    ]
   },
@@ -25,9 +25,34 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "f = LinRange(0.006, 1.0, 100)\n",
+    "num = zeros(100)\n",
     "g0 = g1[1] + g2[1]\n",
-    "T = 100\n",
-    "plotGradient(effects, conc_l, g0, T)"
+    "for (i,v) in enumerate(f)\n",
+    "    ppp = effects[:, 1]\n",
+    "    ppp[4] = v\n",
+    "    num[i] = numcells(ppp, g0, 50)\n",
+    "end\n",
+    "using Plots;\n",
+    "plot(f,num)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "concentration = LinRange(minimum(conc_l), maximum(conc_l), 20)\n",
+    "conc = Vector{Float64}(concentration)\n",
+    "effects2 = getODEparams(ps, conc)\n",
+    "\n",
+    "g0 = g1[1] + g2[1]\n",
+    "xx = LinRange(90.0, 200.0, 200)\n",
+    "@gif for (i,v) in enumerate(xx)\n",
+    "    plotGradient(effects2, concentration, g0, Int(floor(v)), string(i))\n",
+    "    ylims!((-500.0, 500.0))\n",
+    "end every 1\n"
    ]
   },
   {

gemcitabine.ipynb

@@ -9,12 +9,13 @@
     "import Pkg;\n",
     "Pkg.instantiate()\n",
     "using DrugResponseModel\n",
+    "using Plots;\n",
     "\n",
     "# import data from the path\n",
     "conc_l, pop, g2, g1 = setup_data(\"gemcitabine\")\n",
     "\n",
-    "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 4E4)\n",
-    "ps = [10.8326, 1.00369, 1.03549, 2.19808, 0.444854, 1.17305, 0.00923942, 0.122813, 0.483379, 20.7231, 15.4483, 6.60142, 4.18506]\n",
+    "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 6E4)\n",
+    "ps = [26.5608, 0.753368, 1.00515, 1.0, 0.310569, 1.0, 0.00271854, 0.1414, 0.513829, 15.6571, 10.4816, 45.2712, 2.48862]\n",
     "effects = getODEparams(ps, conc_l)"
    ]
   },
@@ -24,9 +25,34 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "f = LinRange(0.0, 0.005, 100)\n",
+    "num = zeros(100)\n",
+    "g0 = g1[1] + g2[1]\n",
+    "for (i,v) in enumerate(f)\n",
+    "    ppp = [0.7533676235831305, 0.310568939766208, 0.0, 0.0, 0.5138291965625418, 15.0, 10.0, 45.0, 2.0]\n",
+    "    ppp[4] = v\n",
+    "    num[i] = numcells(ppp, g0, 150)\n",
+    "end\n",
+    "using Plots;\n",
+    "plot(f,num)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "concentration = LinRange(minimum(conc_l), maximum(conc_l), 20)\n",
+    "conc = Vector{Float64}(concentration)\n",
+    "effects2 = getODEparams(ps, conc)\n",
+    "\n",
     "g0 = g1[1] + g2[1]\n",
-    "T = 100\n",
-    "plotGradient(effects, conc_l, g0, T)"
+    "xx = LinRange(1.0, 150, 100)\n",
+    "@gif for (i,v) in enumerate(xx)\n",
+    "    plotGradient(effects, conc_l, g0, Int(floor(v)), string(i))\n",
+    "    ylims!((-500.0, 500.0))\n",
+    "end every 1"
    ]
   },
   {

lapatinib.ipynb

@@ -9,12 +9,14 @@
     "import Pkg;\n",
     "Pkg.instantiate()\n",
     "using DrugResponseModel\n",
+    "using Plots;\n",
     "\n",
     "# import data from the path\n",
     "conc_l, pop, g2, g1 = setup_data(\"lapatinib\")\n",
     "\n",
-    "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 2E4)\n",
-    "ps = [68.856, 2.17748, 0.0233308, 0.942434, 1.5169, 2.98134, 0.0181415, 0.0359511, 0.531034, 48.6838, 52.3097, 1.66619, 4.44761]\n",
+    "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 5E4)\n",
+    "# ps = [68.856, 2.17748, 0.0233308, 0.942434, 1.5169, 2.98134, 0.0181415, 0.0359511, 0.531034, 48.6838, 52.3097, 1.66619, 4.44761]\n",
+    "ps = [60.0257, 0.925401, 6.87587e-9, 0.792481, 0.368774, 0.750335, 0.00827934, 0.0204481, 0.536002, 20.7355, 13.4557, 0.530622, 2.83636]\n",
     "effects = getODEparams(ps, conc_l)"
    ]
   },
@@ -24,9 +26,34 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "f = LinRange(0.0, 0.82, 100)\n",
+    "num = zeros(100)\n",
+    "g0 = g1[1] + g2[1]\n",
+    "for (i,v) in enumerate(f)\n",
+    "    ppp = effects[:, 1]\n",
+    "    ppp[3] = v\n",
+    "    num[i] = numcells(ppp, g0, 100)\n",
+    "end\n",
+    "using Plots;\n",
+    "plot(f,num)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "concentration = LinRange(minimum(conc_l), maximum(conc_l), 20)\n",
+    "conc = Vector{Float64}(concentration)\n",
+    "effects2 = getODEparams(ps, conc)\n",
+    "\n",
     "g0 = g1[1] + g2[1]\n",
-    "T = 100\n",
-    "plotGradient(effects, conc_l, g0, T)"
+    "xx = LinRange(1.0, 200.0, 200)\n",
+    "@gif for (i,v) in enumerate(xx)\n",
+    "    plotGradient(effects, conc_l, g0, Int(v), string(i))\n",
+    "    ylims!((-200.0, 200.0))\n",
+    "end every 1"
    ]
   },
   {

paclitaxel.ipynb

@@ -9,13 +9,32 @@
     "import Pkg;\n",
     "Pkg.instantiate()\n",
     "using DrugResponseModel\n",
+    "using Plots;\n",
     "\n",
     "# import data from the path\n",
     "conc_l, pop, g2, g1 = setup_data(\"paclitaxel\")\n",
     "\n",
-    "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 4E4)\n",
+    "# costt, ps = optimize_hill(conc_l, g1, g2, maxstep = 6E4)\n",
     "ps = [3.62825, 2.43233, 0.521699, 2.58588, 2.34811, 0.818595, 0.0894338, 0.0482111, 0.469062, 47.4889, 69.1083, 3.90522, 1.12117]\n",
-    "effects = getODEparams(ps, conc_l)"
+    "effects = getODEparams(ps, conc_l);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f = LinRange(0.0, 0.82, 100)\n",
+    "num = zeros(100)\n",
+    "g0 = g1[1] + g2[1]\n",
+    "for (i,v) in enumerate(f)\n",
+    "    ppp = effects[:, 1]\n",
+    "    ppp[4] = v\n",
+    "    num[i] = numcells(ppp, g0, 150)\n",
+    "end\n",
+    "using Plots;\n",
+    "plot(f,num)"
    ]
   },
   {
@@ -24,9 +43,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "concentration = LinRange(minimum(conc_l), maximum(conc_l), 20)\n",
+    "conc = Vector{Float64}(concentration)\n",
+    "effects2 = getODEparams(ps, conc)\n",
+    "\n",
     "g0 = g1[1] + g2[1]\n",
-    "T = 100\n",
-    "plotGradient(effects, conc_l, g0, T)"
+    "xx = LinRange(1.0, 150, 150)\n",
+    "@gif for (i,v) in enumerate(xx)\n",
+    "    plotGradient(effects, conc_l, g0, Int(floor(v)), string(i))\n",
+    "end every 1"
    ]
   },
   {

src/DrugResponseModel.jl

@@ -8,6 +8,7 @@ using Measures
 using BlackBoxOptim
 using LinearAlgebra
 using Base.Threads
+import ExponentialUtilities
 import Calculus
 
 include("importData.jl")
@@ -29,6 +30,9 @@ export setup_data,
     plotUnitSensitivity,
     allSensitivity,
     ODEplot_allPerc,
-    plotGradient
+    plotGradient,
+    predict,
+    numcells,
+    diffCell
 
 end # module

src/Hill.jl

@@ -28,7 +28,6 @@ function residHill(hillParams::Vector, concentrations::Vector, g1::Matrix, g2::M
     return res[]
 end
 
-
 """ Gradient of the cost. """
 function residHillG(hillParams::Vector, concentrations::Vector, g1::Matrix, g2::Matrix)
     # Calculate the continuous parameters with central differencing.
@@ -38,13 +37,12 @@ function residHillG(hillParams::Vector, concentrations::Vector, g1::Matrix, g2::
     return Calculus.finite_difference(hillCost, hillParams)
 end
 
-
 """ Hill optimization function. """
 function optimize_hill(conc_l::Vector, g1::Matrix, g2::Matrix; maxstep = 1E5)
     hillCost(hillParams) = residHill(hillParams, conc_l, g1, g2)
 
-    low = [minimum(conc_l), 1e-9, 1e-9, 0.1, 1e-9, 1e-9, 0.0, 0.0, 0.45, 2, 10, 0, 0]
-    high = [maximum(conc_l), 3.0, 3.0, 10.0, 3.0, 3.0, 1.0, 1.0, 0.55, 60, 180, 50, 50]
+    low = [minimum(conc_l), 1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 0.0, 0.0, 0.35, 2, 10, 0, 0]
+    high = [maximum(conc_l), 1.0, 10.0, 1.0, 10.0, 1.0, 1.0, 1.0, 0.65, 60, 180, 50, 50]
 
     results_ode = bboptimize(
         hillCost;
@@ -60,7 +58,7 @@ end
 
 """ A function to convert the estimated hill parameters back to ODE parameters. """
 function getODEparams(p::Vector, concentrations::Vector{Float64})
-    effects = Matrix{eltype(p)}(undef, 9, 8)
+    effects = Matrix{eltype(p)}(undef, 9, length(concentrations))
 
     # Scaled drug effect
     xx = 1.0 ./ (1.0 .+ (p[1] ./ concentrations) .^ p[4])
@@ -93,7 +91,7 @@ end
 """ Calculate the sensitivity to all parameters. """
 function allSensitivity(params::Vector, conc_l::Vector, g1::Matrix, g2::Matrix)
     b = copy(params)
-    convRange = 10 .^ (range(-1, stop = 1, length = 101))
+    convRange = 10 .^ (range(-0.5, stop = 0.5, length = 101))
     results = zeros(length(convRange), 11)
     paramRanges = zeros(length(convRange), 11)
 
@@ -105,7 +103,7 @@ function allSensitivity(params::Vector, conc_l::Vector, g1::Matrix, g2::Matrix)
     return results, paramRanges
 end
 
-""" Plots the sensitivity for a parameter with a vertical line of the real value of the parameter."""
+""" Plots the sensitivity for a parameter with a vertical line of the real value of the parameter. """
 function plotUnitSensitivity(paramRange, result, realParam, i)
     label = [
         "EC50",
@@ -135,35 +133,34 @@ function plotUnitSensitivity(paramRange, result, realParam, i)
     ylims!((1E2, 1E4))
 end
 
-
 """ Calculate the # of cells in G1 for a set of parameters and T """
-function numcells(params, g0, T)
-    t = LinRange(0.0, 95.5, 192)
-    G1, G2 = predict(params, g0, t, Int(floor(params[6])), Int(floor(params[7])), Int(floor(params[8])), Int(floor(params[9])))
-
-    return G1[T] + G2[T]
+function numcells(params, g0, t::Int)
+    tt = LinRange(0.0, 200, 201)
+    G1, G2 = predict(params, g0, tt, Int(floor(params[6])), Int(floor(params[7])), Int(floor(params[8])), Int(floor(params[9])))
+    @assert(all(x->x>=0, G1), "negative cell number! G1 number is $G1 with this parameter set: $params")
+    @assert(all(x->x>=0, G2), "negative cell number! G2 number is $G2 with this parameter set: $params")
+    @assert(0 <= t <= 200, "you have chosen a time point beyond the simulation, please select 0 <= t <= 200")
+    return G1[t] + G2[t]
 end
 
-
 """ Calculates the gradient with central difference"""
 function diffCell(params, g0, T)
     diffcells(x) = numcells(x, g0, T)
-
-    return Calculus.finite_difference(diffcells, params) / diffcells(params)
+    difs = Calculus.finite_difference(diffcells, params)
+    return difs
 end
 
-
 """ Plot the gradient vs concentrations """
-function plotGradient(effects, concentration, g0, T)
-    dif = zeros(4, 8)
-    for i = 1:8
+function plotGradient(effects, concentration, g0, T, title)
+    dif = zeros(4, length(concentration))
+    for i = 1:length(concentration)
         dif[:, i] = diffCell(effects[:, i], g0, T)[1:4]
     end
     concentrations = log.(concentration)
-    p1 = plot(concentrations, dif[1, :], lw = 2, label = "alpha", xlabel = "log concentration [nM]", ylabel = "gradient of #cells wrt param")
+    p1 = plot(concentrations, dif[1, :], lw = 2, label = "alpha", xlabel = "log concentration [nM]", ylabel = "gradient of #cells wrt param", title=title)
     plot!(concentrations, dif[2, :], lw = 2, label = "beta")
-    p2 = plot(concentrations, dif[3, :], label = "gamma1", lw = 2, xlabel = "log concentration [nM]", ylabel = "gradient of #cells wrt param")
+    p2 = plot(concentrations, dif[3, :], label = "gamma1", lw = 2, xlabel = "log concentration [nM]", ylabel = "gradient of #cells wrt param", title=title)
     plot!(concentrations, dif[4, :], lw = 2, label = "gamma2")
     plot(p1, p2, layout = (1, 2))
     plot!(size = (800, 400), margin = 0.4cm, dpi = 100)
-end
+end