Skip to content

Instantly share code, notes, and snippets.

@saschatimme

saschatimme/Solving a chemical reaction network with HomotopyContinuation.jl.ipynb

Created January 25, 2019 13:37
Show Gist options
  • Save saschatimme/38fc2bbe641cfd7f3985c0169f6efed6 to your computer and use it in GitHub Desktop.
Save saschatimme/38fc2bbe641cfd7f3985c0169f6efed6 to your computer and use it in GitHub Desktop.
Download ZIP
Solving a chemical reaction network with HomotopyContinuation.jl.ipynb
Display the source blob
Display the rendered blob
Raw
Solving a chemical reaction network with HomotopyContinuation.jl.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** Disclaimer**: This demonstration works best with a soon to be relased version of [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org).\n",
"You can install this with"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import Pkg\n",
"Pkg.add(\"https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl.git#releases/0.5\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Solving a chemical reaction network with HomotopyContinuation.jl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to solve a bio-chemical reaction newtwork.\n",
"The system comes from biochemical reactions ([published here](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005267)). Basically it models how the bacteris B. subtilis produces a stress response protein, sigmaB, in response to some stress input. This input to the system is the concentration of a phosphatase (a kind of protein).\n",
"\n",
"Here is the system with parameters:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10-element Array{Polynomial{true,Float64},1}:\n",
" -x₁²x₇p₁ - x₁²p₁p₂₀ + x₇p₁₈p₁₉p₂₁ - x₁x₇p₁₇ - x₁p₁₇p₂₀ + 2.0x₂x₇p₂ + 2.0x₂p₂p₂₀ + p₁₈p₂₁ \n",
" 0.5x₁²p₁ - x₂x₄p₄ - x₂x₇p₆ - x₂p₂ - x₂p₁₇ + x₃p₉ + x₃p₁₄ + x₈p₁₁ \n",
" x₂x₄p₄ - x₃x₄p₅ - x₃x₇p₁₂ + x₄x₈p₇ - x₃p₉ - x₃p₁₄ - x₃p₁₇ + x₅p₁₀ + x₅p₁₅ \n",
" -x₂x₄x₇p₄ - x₂x₄p₄p₂₀ - x₃x₄x₇p₅ - x₃x₄p₅p₂₀ + x₃x₇²p₁₂ + x₃x₇p₁₂p₂₀ - x₄x₇x₈p₇ - x₄x₈p₇p₂₀ + x₇p₁₈p₁₉p₂₂ + x₃x₇p₉ + x₃p₉p₂₀ - x₄x₇p₁₇ - x₄p₁₇p₂₀ + x₅x₇p₁₀ + x₅p₁₀p₂₀ + x₇x₉p₁₆ + x₉p₁₆p₂₀ + p₁₈p₂₂\n",
" x₃x₄p₅ - x₅p₁₀ - x₅p₁₅ - x₅p₁₇ \n",
" -x₆x₁₀p₈ + x₃p₁₄ + x₅p₁₅ - x₆p₁₇ + x₉p₁₃ \n",
" -x₂x₇²p₆ - x₂x₇p₆p₂₀ - x₃x₇²p₁₂ - x₃x₇p₁₂p₂₀ + x₄x₇x₈p₇ + x₄x₈p₇p₂₀ - x₇²p₁₇ + x₇x₈p₁₁ - x₇p₁₇p₂₀ + x₇p₁₈p₁₉ + x₈p₁₁p₂₀ + p₁₈ \n",
" x₂x₇p₆ + x₃x₇p₁₂ - x₄x₈p₇ - x₈p₁₁ - x₈p₁₇ \n",
" x₆x₁₀p₈ - x₉p₁₃ - x₉p₁₆ - x₉p₁₇ \n",
" x₉ + x₁₀ - p₂₃ "
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using HomotopyContinuation, DynamicPolynomials\n",
"\n",
"@polyvar x[1:10] p[1:23]\n",
"\n",
"f = [\n",
" (-1 * p[17] * x[1] + -2 * p[1] * (x[1] ^ 2 / 2) + 2 * p[2] * x[2])*(p[20] + x[7]) + p[21] * p[18] * ((1 + p[19] * x[7])),\n",
" -1 * p[17] * x[2] + p[1] * (x[1] ^ 2 / 2) + -1 * p[2] * x[2] + -1 * p[4] * x[2] * x[4] + p[9] * x[3] + p[14] * x[3] + -1 * p[6] * x[2] * x[7] + p[11] * x[8],\n",
" -1 * p[17] * x[3] + p[4] * x[2] * x[4] + -1 * p[9] * x[3] + -1 * p[5] * x[3] * x[4] + p[10] * x[5] + -1 * p[14] * x[3] + p[15] * x[5] + p[7] * x[8] * x[4] + -1 * p[12] * x[3] * x[7],\n",
" (-1 * p[17] * x[4] + -1 * p[4] * x[2] * x[4] + p[9] * x[3] + -1 * p[5] * x[3] * x[4] + p[10] * x[5] + -1 * p[7] * x[8] * x[4] + p[12] * x[3] * x[7] + p[16] * x[9])*(p[20] + x[7]) + p[22] * p[18] * (1 + p[19] * x[7]),\n",
" -1 * p[17] * x[5] + p[5] * x[3] * x[4] + -1 * p[10] * x[5] + -1 * p[15] * x[5],\n",
" -1 * p[17] * x[6] + p[14] * x[3] + p[15] * x[5] + -1 * p[8] * x[6] * x[10] + p[13] * x[9],\n",
" (-1 * p[17] * x[7] + -1 * p[6] * x[2] * x[7] + p[11] * x[8] + p[7] * x[8] * x[4] + -1 * p[12] * x[3] * x[7])*(p[20] + x[7]) + p[18] * (1 + p[19] * x[7]),\n",
" -1 * p[17] * x[8] + p[6] * x[2] * x[7] + -1 * p[11] * x[8] + -1 * p[7] * x[8] * x[4] + p[12] * x[3] * x[7],\n",
" -1 * p[17] * x[9] + p[8] * x[6] * x[10] + -1 * p[13] * x[9] + -1 * p[16] * x[9],\n",
" x[10] + x[9] - p[23]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are interested in real solutions where each entry is non-negative and we want to solve this system for many different choices of parameters."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By Bezout's theorem we know that the polynomial system $f$ has at most `prod(maxdegree,f) == 1728` solutions.\n",
"But the root count is probably much less.\n",
"By using algebraic geometry we know that for almost all choices of $p \\in \\mathbb{C}^{23}$ the system $f(x;p)$ has same number of solutions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A strategy to solve many instance of our problem is therefore to first solve a system for random `p`. Then we can track the solutions of the system from one parameter to another one."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"AffineResult with 1728 tracked paths\n",
"==================================\n",
"• 44 non-singular finite solutions (0 real)\n",
"• 0 singular finite solutions (0 real)\n",
"• 1530 solutions at infinity\n",
"• 154 failed paths\n",
"• random seed: 413742\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Plug in some random parameters and solve them\n",
"p_rand = randn(ComplexF64, 23)\n",
"f_rand = subs.(f, Ref(p=> p_rand))\n",
"\n",
"result_rand = solve(f_rand, report_progress=false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the system has generically only **44** solutions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now can use this to find solutions to some specific set of parameters (which also are physically meaningful):"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"AffineResult with 44 tracked paths\n",
"==================================\n",
"• 42 non-singular finite solutions (11 real)\n",
"• 2 singular finite solutions (0 real)\n",
"• 0 solutions at infinity\n",
"• 0 failed paths\n",
"• random seed: 819499\n"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# parameters we are interested in\n",
"p₀ = [3600, 18, 18, 3600, 3600, 3600, 1800, 3600, 18, 18, 18, 1800, 18, 36, 11, 180, 0.7, 0.4, 30, 0.2, 4, 4.5, 0.4]\n",
"\n",
"# Now make a parameter homotopy from our generic start system\n",
"result_p₀ = solve(f, solutions(result_rand), parameters=p, p₁=p_rand, p₀=p₀)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we only want positive real solutions, lets make a helper function:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"positive_real_solutions (generic function with 1 method)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function positive_real_solutions(result)\n",
" filter(realsolutions(result)) do x\n",
" all(xᵢ -> xᵢ ≥ 0, x)\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"No we can find the solution:"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10-element Array{Float64,1}:\n",
" 0.11584484227605159 \n",
" 0.9455229982635467 \n",
" 1.4935460651014791 \n",
" 0.014733757161745744\n",
" 2.6673387918892226 \n",
" 15.827943586937968 \n",
" 0.03777228097321519 \n",
" 5.088785714445814 \n",
" 0.3986099864060771 \n",
" 0.001390013593923053"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"positive_real_solutions(result_p₀)[1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or some slightly different parameters:"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1-element Array{Array{Float64,1},1}:\n",
" [0.10192, 0.43226, 1.78391, 0.0261887, 5.66282, 21.6704, 0.106577, 7.71679, 0.398984, 0.00101621]"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"q₀ = copy(p₀)\n",
"q₀[14] = 14.0\n",
"\n",
"result_q₀ = solve(f, solutions(result_rand), parameters=p, p₁=p_rand, p₀=q₀)\n",
"positive_real_solutions(result_q₀)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.3",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment