Numerical integration methods for large-scale biophysical simulations

Simulations of biophysical systems inevitably include steps that correspond to time integrations of ordinary differential equations. These equations are often related to enzyme action in the synthesis and destruction of molecular species, and in the regulation of transport of molecules into and out of the cell or cellular compartments. Enzyme action is almost invariably modeled with the quasi-steady- state Michaelis–Menten formula or its close relative, the Hill formula: this description leads to systems of equations that may be stiff and hard to integrate, and poses unusual computational challenges in simulations where a smooth evolution is interrupted by the discrete events that mark the cells’ lives, like initiation and termination of DNA synthesis or mitosis. These discrete events lead to abrupt parameter changes and to variable system size. This is the case of a numerical model (Virtual Biophysics Lab – VBL) that we are developing to simulate the growth of three-dimensional tumor cell aggregates (spheroids). The underlying cellular events have characteristic timescales that span approximately 12 orders of magnitude, and thus the program must be robust and stable. Moreover the program must be able to manage a very large number of equations (of the order of 107–108), and finally it must be able to accept frequent modifications of the underlying theoretical model. Here we study the applicability of known integration methods to this context – quite unusual from the point of view of the standard theory of differential equations, but extremely relevant to biophysics – and we describe the results of numerical tests in situations similar to those found in actual simulations.