Solving nonlinear coupled reaction diffusion equations (reaction networks) in FiPy

[Alokendra]

Hi,

I am a new potential FiPy user. I have been looking for a python library
that can solve a system of coupled reaction-diffusion equations used in
models of biological pathways. I have been using a software called VCell (
http://vcell.org/) which can solve a wide range of problems such as
deterministic and stochastic system of PDEs and use a standard file
format called SBML (http://sbml.org/Basic_Introduction_to_SBML). But the
VCell interface is quite constrained there and it does not allow easy
parallelization so I was looking for an alternative.

I was wondering how you go about solving a such system of equations with
drift component as shown below in FiPy. This for one particular species
A of the system and I have similar equations for other species like B
and C. The concentration of the species satisfy some common mass-
conservation relation. The equations are mostly coupled through the
source/reaction term which can be nonlinear Hill-type functions.

Equation:
$$\frac{\partial C_A}{\partial t} = D_A\nabla^2C_A - V_A\nabla C_A -
(k_fC_AC_B-k_rC_C)$$

I have read the tutorial and I think I understand how to setup the
diffusion and convection terms but I am not so sure about the reaction
terms which can be coupled nonlinear functions. Do I need to use some
variation of `ImplicitSource` term?

Also what is the most efficient way to get Fipy read information such as
rate constants, diffusion coefficients and drift diffusion coefficients
for many such species? Because of the nonlinearity in the reaction term
it cannot be vectorized in straightforward manner.

Regards,
Alok

Alok Ghosh
Graduate Student
University of Pennsylvania

[Guyer, Jonathan E. Dr. (Fed)]

Post by Alokendra
I am not so sure about the reaction
terms which can be coupled nonlinear functions. Do I need to use some
variation of `ImplicitSource` term?

Yes, these terms are appropriate for ImplicitSourceTerms, e.g., I would probably write $(k_fC_AC_B-k_rC_C)$ as

ImplicitSourceTerm(coeff=k_f * C_A, var=C_B) - ImplicitSourceTerm(coeff=k_r, var=C_C)

The first term could also be written

ImplicitSourceTerm(coeff=k_f * C_B, var=C_A)

but that would result in less coupling between the equations, which would generally slow convergence. The first form might not be stable, though, so I would have to experiment if I was doing this.

Post by Alokendra
Also what is the most efficient way to get Fipy read information such as
rate constants, diffusion coefficients and drift diffusion coefficients
for many such species? Because of the nonlinearity in the reaction term
it cannot be vectorized in straightforward manner.

That depends on the form the coefficients take wherever you're getting them from.

I would argue that one of the biggest advantages of FiPy over some other PDE solvers is that you define your problem with the Turing-complete Python programming language and not some dedicated scripting language or batch input file. You can do anything that Python can do to read and parse files.