# Solving Catalyst/Enzyme Diffusion Equation

Evening all,

I am attempting to solve a second order differential equation to get an equation of the substrate concentration as a function of r (i.e. S(r)). I have developed this equation, $\frac{d^2S}{dr^2}+\frac{2}{r}\frac{dS}{dr}-\frac{V_mS}{K+S}*\frac{1}{D}=0$, where $V_m, D, K$ are known values to me based on research and I also have two boundary conditions as it is a symmetrical spherical enzyme. Thus, $\frac{dS}{dr} = 0, r=0$ and $S=0.00001, R=0.001$ I am very new to Mathematica and have been attempting to solve this equation in the following manner:

DSolve[{S''[r] + (2/r)*S'[r] - (V*S[r])/(K + S[r])*1/D == 0, S' == 0,
S[0.001] == 1.10*10^-4}, S[r], r]


I have not been successful in getting a solution as S(r) and am unsure where to begin on fixing the issue. Any help is greatly appreciated!

• Could use NDSolve with explicit values for the parameters. But it seems you will need to move the initial point to something positive as it otherwise claimes a singularity at the origin. – Daniel Lichtblau Sep 17 '17 at 21:43
• Sorry, completely new to the software. How would I go about using NDSolve for this equation? – JuliusDariusBelosarius Sep 17 '17 at 22:47
• What I had in mind was exactly covered in the response by @zhk. – Daniel Lichtblau Sep 18 '17 at 16:35

## 2 Answers

DSolve is unable to find an exact solution, so we go for numerical one using NDSolve.

You are not suppose to use built-in functionality of Mathematica as a parameter.

I have changed D, which is used for derivative to D1. Apart from that, for r=0, you are facing 1/0. So I choose the starting value to be 10^-4.

V1 = 1; K1 = 1; D1 = 1;

Eq = S''[r] + (2/r)*S'[r] - (V1*S[r])/(K1 + S[r])*1/D1 == 0

sol = NDSolve[{Eq, S'[10^-4] == 0, S[0.001] == 1.10*10^-4}, S[r], r]

Plot[S[r] /. sol, {r, 10^-4, 1}]


You can also use ParametricNDSolve, e.g.:

pn = ParametricNDSolve[{s''[r] + 2 s'[r]/r - v  s[r]/(d (k + s[r])) ==
0, s[0.001] == 1.10 10^(-4), s'[0.000001] == 0},
s, {r, 0.0001, 1}, {v, d, k}]
Manipulate[
Plot[Evaluate[s[v, d, k][r]] /. pn, {r, 0.001, 1},
PlotRange -> {0, 0.0004}], {v, 1, 10}, {d, 1, 10}, {k, 1, 10}] 