Hybrid PDE-Agent Models for Cell-Tissue Systems

Why Couple Fields and Cell Agents?

In cell–tissue systems, different parts of the biology call for different representations. An extracellular mediator may diffuse through tissue and have a concentration at every spatial location, while each cell may need its own state for receptor abundance, internalized material, or phenotype. A hybrid PDE-agent model represents the extracellular substance as a continuous field and each cell as a discrete agent, then specifies how material passes between them.

The coupling is not just a matter of reading a field value at each cell. The cell’s uptake or release must also appear in the field equation with the opposite sign and a spatial distribution that transfers the same amount. The definitions below make this exchange explicit. The model is a general construction, not a claim that one particular transport law applies to every cell type or substance.

Definition: Let \(c(\mathbf{x},t)\) be extracellular concentration, in amount per tissue volume, at location \(\mathbf{x}\) and time \(t\). For cell \(i\), let \(m_i(t)\) be the amount of the modeled substance inside that cell, and let \(z_i(t)=m_i(t)/V_i\) be its intracellular concentration for fixed cell volume \(V_i\). Define \(J_i(t)\) as the net rate of transfer from the extracellular space into cell \(i\), in amount per time. Thus \(J_i>0\) means uptake by the cell, and \(J_i<0\) means release into the extracellular space.

Specify the Exchange and Its Spatial Footprint

A simple membrane law makes exchange proportional to the difference between the local extracellular concentration and the intracellular concentration adjusted for partitioning. Let \(P_i\) be membrane permeability (length per time), \(A_i\) the exchange area (area), and \(K_i\) a dimensionless partition ratio. The law is

J_i=P_iA_i\left(\bar{c}_i-\frac{z_i}{K_i}\right), \qquad \bar{c}_i(t)=\int_{\Omega}W_i(\mathbf{x},t)c(\mathbf{x},t)\,d\mathbf{x}.

Here \(\Omega\) is the tissue domain and \(W_i\) is a nonnegative spatial weighting function centered on or near cell \(i\). It is normalized so that \(\int_{\Omega}W_i(\mathbf{x},t)\,d\mathbf{x}=1\). Consequently, \(\bar{c}_i\) is a local weighted extracellular concentration, not an unweighted tissue average. The units of \(W_i\) are inverse volume. For this linear law, zero flux occurs when \(z_i=K_i\bar{c}_i\). A different biological transport mechanism may require a different flux law, but its units and sign convention must remain explicit.

The equal-and-opposite field contribution is \(-J_iW_i\). Because \(J_i\) has units of amount per time and \(W_i\) has units of inverse volume, this term has units of concentration per time, as required in a concentration PDE. One common coupled model is

\frac{\partial c}{\partial t} = D\nabla^2 c+R(\mathbf{x},t)-k c -\sum_{i=1}^{N}J_i(t)W_i(\mathbf{x},t), \qquad \frac{dm_i}{dt}=J_i(t)-v_i(t).

In this expression, \(D\) is the extracellular diffusion coefficient, \(R\) is any external production rate per tissue volume, \(k\) is a first-order extracellular loss rate, and \(v_i\) is intracellular removal or transformation in amount per time. If \(v_i\) represents degradation inside a cell, it removes the modeled substance from the combined extracellular-plus-intracellular inventory. It is not an exchange flux and must not be added as an extracellular sink unless the biology specifically returns a product represented by the field.

A localized kernel avoids placing the entire cell exchange at an ideal point, which can create an unphysical singularity in a continuum field. The support and shape of \(W_i\) should correspond to the spatial resolution and intended biological interpretation. In a spatial discretization, the kernel weights must be mapped to grid cells so their volume-weighted sum is one; otherwise the grid receives too much or too little material.

Worked Example: Instantaneous Uptake Rate

Hypothetical values. A cell has \(P_i=2.0\times10^{-7}\,\mathrm{m/s}\), \(A_i=1.0\times10^{-9}\,\mathrm{m^2}\), \(K_i=1\), local extracellular concentration \(\bar{c}_i=4.0\,\mathrm{\mu mol/L}\), and intracellular concentration \(z_i=1.0\,\mathrm{\mu mol/L}\). Calculate the instantaneous membrane flux.

First convert concentration because \(P_iA_i\) is in \(\mathrm{m^3/s}\): \(1\,\mathrm{\mu mol/L}=1.0\times10^{-3}\,\mathrm{mol/m^3}\). Therefore \(\bar{c}_i-z_i/K_i=3.0\times10^{-3}\,\mathrm{mol/m^3}\), and \(P_iA_i=2.0\times10^{-16}\,\mathrm{m^3/s}\).

J_i=(2.0\times10^{-16}\,\mathrm{m^3/s}) (3.0\times10^{-3}\,\mathrm{mol/m^3}) =6.0\times10^{-19}\,\mathrm{mol/s}.

The positive sign indicates uptake at the instant represented by these concentrations. This calculation gives a rate at that state. It does not, by itself, establish how much material will enter the cell over a later finite interval: as exchange proceeds, the intracellular and local extracellular concentrations may change, changing \(J_i\).

Mass-Balance Justification

The opposite signs in the two equations are required by conservation during exchange. To see this, take a fixed tissue domain with no-flux boundaries, so diffusion contributes no net amount across its boundary. Integrate the field equation over the domain. Normalization gives \(\int_{\Omega}W_i\,d\mathbf{x}=1\), so each cell’s contribution to the rate of change of extracellular amount is \(-J_i\). Meanwhile, the rate of change of intracellular amount for that cell is \(J_i-v_i\). Adding the extracellular and intracellular rates cancels exchange exactly.

\frac{d}{dt}\left(\int_{\Omega}c\,d\mathbf{x}+\sum_{i=1}^{N}m_i\right) = \int_{\Omega}R\,d\mathbf{x} -k\int_{\Omega}c\,d\mathbf{x} -\sum_{i=1}^{N}v_i.

This balance assumes that \(c\) is amount per volume, \(m_i\) is amount, each \(W_i\) is normalized over the modeled domain, and the boundary condition contributes no net diffusive flux. With a boundary that permits transport, its net flux must be included in the balance. The equation also distinguishes exchange from production and degradation: exchange redistributes the modeled substance between compartments, while \(R\), \(k c\), and \(v_i\) change its combined inventory.

Key result: A consistent exchange term uses the same \(J_i\) in both model components: \(+J_i\) in the cell’s amount balance and \(-J_iW_i\) in the field equation. A normalized \(W_i\) ensures that the field loses exactly \(J_i\) amount per time when the cell takes up material.

Worked Example: Finite-Time Exchange in a Closed Two-Compartment System

Hypothetical setup. Consider one well-mixed extracellular compartment of volume \(V_e=1.0\times10^{-12}\,\mathrm{m^3}\) coupled to one cell of volume \(V_i=1.0\times10^{-15}\,\mathrm{m^3}\). There is no diffusion, production, degradation, or other loss. Let \(J=g(c-z)\), with \(g=P_iA_i=1.0\times10^{-15}\,\mathrm{m^3/s}\), \(K_i=1\), initial extracellular concentration \(c(0)=4.0\times10^{-3}\,\mathrm{mol/m^3}\), and initial intracellular concentration \(z(0)=1.0\times10^{-3}\,\mathrm{mol/m^3}\). Find the amount transferred into the cell in \(0.100\,\mathrm{s}\).

The coupled equations are \(V_e\,dc/dt=-J\) and \(V_i\,dz/dt=J\). Total amount \(V_ec+V_iz\) is constant. The concentration difference \(d=c-z\) obeys

\frac{dd}{dt} = -g\left(\frac{1}{V_e}+\frac{1}{V_i}\right)d, \qquad d(t)=d(0)e^{-\lambda t}, \qquad \lambda=g\left(\frac{1}{V_e}+\frac{1}{V_i}\right).

Substitution gives \(\lambda=(1.0\times10^{-15})(1.0\times10^{12}+1.0\times10^{15})\,\mathrm{s^{-1}}=1.001\,\mathrm{s^{-1}}\). Integrating \(J=g d(t)\) over the interval gives the transferred amount:

\Delta m_i = \int_0^{0.100}J(t)\,dt = \frac{g\,d(0)}{\lambda}\left(1-e^{-\lambda(0.100\,\mathrm{s})}\right) \approx 2.85\times10^{-19}\,\mathrm{mol}.

For the arithmetic, \(g\,d(0)/\lambda=(1.0\times10^{-15})(3.0\times10^{-3})/1.001\approx2.997\times10^{-18}\,\mathrm{mol}\), and \(1-e^{-0.1001}\approx0.09525\). Thus the cell gains approximately \(2.85\times10^{-19}\,\mathrm{mol}\), while the extracellular compartment loses the same amount. Simply multiplying the initial rate \(J(0)=3.0\times10^{-18}\,\mathrm{mol/s}\) by \(0.100\,\mathrm{s}\) gives \(3.0\times10^{-19}\,\mathrm{mol}\), an explicit constant-rate estimate, not the exact finite-time transfer. The exact result is smaller because the concentration difference and flux decline during exchange.

Mapping Exchange onto a Spatial Grid

In computation, a cell’s flux must be assigned to field locations without changing its integrated amount. Suppose a cell’s normalized kernel assigns weights \(w_{ij}\) to grid volumes \(V_j\), where \(w_{ij}\geq0\) and \(\sum_jw_{ij}=1\). The field concentration rate in grid volume \(j\) receives \(-J_iw_{ij}/V_j\). Multiplying that rate by \(V_j\) recovers the amount rate deposited in that volume.

Worked Example: Checking a Discrete Exchange Map

Hypothetical values. A cell has instantaneous uptake \(J_i=4.0\times10^{-18}\,\mathrm{mol/s}\). Its kernel allocates 25% of the exchange to voxel 1 and 75% to voxel 2. Each voxel has volume \(2.0\times10^{-12}\,\mathrm{m^3}\). Calculate the field concentration-rate contributions and verify the total loss.

For voxel 1, the allocated amount rate is \(J_iw_{i1}=(4.0\times10^{-18})(0.25)=1.0\times10^{-18}\,\mathrm{mol/s}\). Its concentration-rate contribution is therefore \(-1.0\times10^{-18}/(2.0\times10^{-12})=-5.0\times10^{-7}\,\mathrm{mol\,m^{-3}\,s^{-1}}\).

For voxel 2, the allocated amount rate is \((4.0\times10^{-18})(0.75)=3.0\times10^{-18}\,\mathrm{mol/s}\). The concentration-rate contribution is \(-3.0\times10^{-18}/(2.0\times10^{-12})=-1.5\times10^{-6}\,\mathrm{mol\,m^{-3}\,s^{-1}}\).

Multiplying each concentration rate by its voxel volume gives field losses of \(1.0\times10^{-18}\) and \(3.0\times10^{-18}\,\mathrm{mol/s}\), respectively. Their sum is \(4.0\times10^{-18}\,\mathrm{mol/s}\), exactly the cell’s uptake rate. The negative signs in the field equation indicate decreases in extracellular concentration.

Assumptions, Diagnostics, and Alternatives

The linear membrane law assumes that a permeability-area product and an effective partition ratio adequately describe net exchange over the modeled range. It also assumes that the concentrations used in the law are well-defined at the cell–field interface or are represented acceptably by the chosen local average. Fixed \(V_i\) is needed when converting between cellular amount and concentration by \(z_i=m_i/V_i\); if cell volume changes, the amount balance should be written first, and concentration changes should account for volume changes. The tissue field is assumed to be meaningful at the spatial and temporal resolution used.

Useful checks include the following:

Units. Verify that every field term has concentration per time units, every agent amount equation has amount per time units, and the flux used in both equations has amount per time units.
Kernel normalization. Check the integral normalization in a continuous model or the volume-weighted grid normalization in a discretization. A sum of weights alone is not sufficient if the weights represent densities rather than dimensionless allocations.
Sign and conservation. For a closed, exchange-only test, positive cell uptake must produce an equal extracellular loss. A changing total amount signals a sign, normalization, or implementation error.
Time-step behavior. Large explicit time steps can overshoot, drive concentrations negative, or transfer more material than is available. Compare time-step refinements and, for stiff exchange, consider an implicit or analytically integrated exchange update. Do not assume that the initial flux stays constant across a finite step unless deliberately making and labeling that approximation.
Spatial resolution. If the exchange kernel is narrower than the grid, its mapped profile may depend strongly on grid placement. Check whether the local concentration sampled by the cell and the resulting uptake are stable under reasonable grid refinement.

A continuum-only model can be appropriate when individual cell differences and discrete events are not needed; it replaces cell-specific states with averaged cellular quantities. An agent-only model can represent individual cells but may not naturally resolve smooth extracellular gradients without additional spatial rules. A hybrid model retains individual cellular states while solving a continuum transport equation, at the cost of defining and verifying their exchange. It is not automatically more realistic: if the data cannot identify cell-level behavior or spatial heterogeneity, its additional structure may be unsupported.

Common pitfall: A modeler computes the field loss as \(-J_i\) at a grid point and the cell gain as \(+J_i\), then treats the field term as a concentration rate. This is dimensionally inconsistent: \(J_i\) is amount per time, whereas a concentration equation needs amount per volume per time. Distribute the flux with a normalized kernel, or divide an allocated amount rate by the corresponding spatial volume, and verify that integrating the field loss recovers \(-J_i\).

Check Your Understanding

Answer each question, then compare with the concise solutions below.

If \(J_i<0\), what does the sign mean for the cell and the extracellular field?
Why must the spatial kernel be normalized, and what units does it have in a continuous three-dimensional model?
A model includes \(+J_i\) in the cell equation but \(-2J_iW_i\) in the field equation, with normalized \(W_i\). What happens to exchange-only mass balance?
For a closed system, name two processes in the displayed model that can change the combined extracellular-plus-cellular amount.
Why is multiplying \(J_i(0)\) by a finite time interval generally not an exact calculation of transferred amount?

Solutions. (1) The cell releases material; the field receives it because its exchange term is \(-J_iW_i>0\). (2) Normalization makes the integrated field contribution equal to the intended amount rate. In three dimensions \(W_i\) has units of inverse volume. (3) The field loses twice the amount the cell gains, so the modeled combined amount decreases at rate \(J_i\) for positive \(J_i\); exchange is not conservative. (4) External production \(R\), extracellular degradation \(k c\), and intracellular removal \(v_i\) can change the combined amount; any two suffice. (5) The flux depends on concentrations that generally change during the interval. The product \(J_i(0)\Delta t\) is only a constant-initial-rate estimate, such as an explicit-Euler approximation, unless the flux is in fact constant.