Mathematical model of the post-ablation enhancement zone as a tissue-level oedematic response

Abstract

Purpose: A hyperdense rim is commonly observed at the periphery of ablation zones during post-ablation imaging (e.g. ultrasound) in tumours. A mathematical model has been developed here to investigate the occurrence of this enhanced rim, caused by the ablated cells, giving an indication of the location of the final ablation region.

Materials and methods: The enhanced rim has been assumed here to be due to a tissue-level oedematic response of viable cells, which necessitated coupling multiple modelling elements in a spatially distributed system: thermal cell death, tissue-state dependent ion concentration dynamics, ion transport in the extracellular space, and osmotic cell volume regulation.

Results: In response to the imposed temperature function, an ablation zone was predicted, distinguishing the tissue state between ‘dead’ and ‘alive’. A disturbance in intracellular/extracellular ion concentrations was induced due to ion redistribution, which acted as an osmotic stress and contributed to significant cell swelling in a thin rim at the periphery of the ablation zone. It was also found that the rim size only changed slightly with varying lesion size, in response to different temperature profiles.

Conclusions: The study presents a novel mathematical model to understand the enhanced rim surrounding the ablation zone by assuming tissue-level cell oedema as the primary potential cause. The model links the direct response to thermal injury to an observable secondary response, which could be of clinical value in that the location of this bright ring could potentially be used for more accurate determination of the extent of the ablation zone.

Keywords:

• Cellular oedema
• enhancement rim
• mathematical model
• physiological response
• thermal ablation

Introduction

Minimally invasive hyperthermic ablations with guidance from imaging play an increasing role in cancer treatment. Radiofrequency ablation (RFA) and microwave ablation (MWA) are two of the most commonly used treatments. The underlying mechanism of thermal therapies is that irreversible cell damage is induced in response to hyperthermic temperatures (≥50 °C), with cells ultimately undergoing coagulative necrosis [1]. Incomplete ablation is one major cause of failure of thermal therapies, and this can often be attributed to the heat-sink effect of the tumour vasculature [1] and the lack of experience of the radiotherapists [2]. Mathematical modelling provides a valuable tool to assist in treatment planning, with considerable effort being spent on the development of bio-heat transfer models [3] and on modelling thermal damage processes and cell death [4].

Clinically, the response of thermal treatments is assessed by post-treatment imaging, often contrast-enhanced computed tomography (CT), with the appearance and size of ablation regions being of primary concern. In addition, a rim of enhancement has been observed at the periphery of the ablation region in different phases, ranging from immediately after treatment to several months, indicating an increase in density of tissue [5–7]. It has been postulated to be a physiological response to thermal injury, which initially results from hyperaemia, and subsequently from fibrosis and giant cell reaction [5]. The hyperintense rim is also a common finding from post-ablative MR-diffusion-weighted imaging, which indicates a significant drop in the apparent diffusion coefficient (ADC), one that might be correlated to oedema as a reaction to thermal stress [8]. Cellular oedema is also hypothesised to contribute to its presence [9]. Analogous to cytotoxic oedema in the field of neuroradiology [10], oedema may be induced in response to the osmotic stress which results from the movement of active ions in the extracellular space due to impairment of trans-membrane ion pumps as a consequence of thermal cell death. In this regard, the presence of the rim is closely tied to the final lesion size, which could be of clinical value in terms of providing supplementary information about and evaluation of the ablation volume.

This study presents a novel mathematical model to capture the occurrence of cellular swelling at the periphery of ablation zones, which was presumed to be the dominant mechanism for the rim of enhancement. As there are few published studies for modelling such behaviour, an initial step was taken to develop a mathematical model that integrates three essential components: (1) a three-state thermal cell death model [11] was adopted to predict the size of the ablation region, (2) an ion-based model has been included to account for the kinetics of ionic trans-membrane transport and regulation of cell volume controlled by osmotic processes, (3) ion transport in the extracellular space plays a key role in the movement of the effects of thermal cell death on trans-membrane ionic kinetics in the lesion zone to its neighbouring region, generating a redistribution of ion concentrations and thus imposing osmolality stresses. By having all the elements assembled in a well-controlled computational modelling set-up, clear explanations can be provided on how cell swelling is induced as the physiological response of cells at the periphery of ablation zones.

Methods

Figure 1 illustrates the components of the model and their interplay contributing to cellular oedema. A mechanical thermal cell death model has been used to predict the tissue state in terms of the spatial distribution of cell viability in response to hyperthermia damage. Up to now, most of the modelling efforts on cellular oedema have been made in the cardiac or nervous systems [12–16]. Despite very sparse literature being available with respect to the tumour environment, a generic Hodgkin–Huxley formulation for ionic currents has been used here to track ionic transport dynamically. Sodium, potassium and chloride were specified due to their well-established transport (pump-leak) modes, and the cell volume control was regulated by transmembrane osmotic pressure [15,16]. Therefore, the ion-based model has been presented in general terms without referring to any specific cell type, the details related to which (e.g. the specific ion species and transport systems involved) could and should be examined individually in future.

Figure 1. Schematic diagram of model of cell oedema post-ablation.

An alteration in the dynamics of ion concentrations is caused by the transition of cell states. When embedded in a spatially distributed system, any disturbance in local behaviour is not isolated; rather, it is transported to its neighbourhood by ion movement. In this way any deviation in ion concentrations from their physiological levels can be potentially presented as a stress signal, which induces an osmotic response, with a potential increase in cell volume being observed.

Modelling dynamics of intracellular ionic concentrations

In the current work, concentrations of ions have been denoted by Na_e/i, [K]_e/i and [Cl]_e/i with subscripts ‘e’ and ‘i’ referring to the extracellular and intracellular spaces respectively. In normothermia, intracellular ion quantities are balanced by ion currents across the cell membrane through ion channels (i_Na, i_K, i_Cl), and the ion pump (the Na-K-ATP pump with exchange of Na⁺ and K⁺ ions in the ratio 3:2, i_NaK), as shown in Figure 1 and described in the following equations, based on the pump-leak mechanism [15]: (1) (2) (3) where F is Faraday’s constant, and v_c is the cell volume (dynamically regulated by osmolality). The ion currents from passive ion transport were assumed to be described by a Hodgkin–Huxley-type model, in terms of membrane conductances, the membrane potential (E_m) and the Nernst potentials (v_ion): (4) (5) (6)

For Na⁺ and K⁺, the membrane conductances were a combination of the background leak conductances (g_Na^l and g_K^l) and voltage-dependent gated conductances as a function of voltage-dependent gating variables taken from the Hodgkin–Huxley model [16]. These were assumed to be at quasi-steady state (m_∞, h_∞ and n_∞) due to the very fast timescales of the gating of ions, as given in Table 1. Though the time constants of channel gating (i.e. ion channel dynamics) were more sensitive to the temperature change [17], the variables at quasi-steady state were of interest. Thus, it has been assumed that the impact of temperature on ion dynamics can be neglected in the current study.

Table 1. Gating variables at quasi-steady state (m_∞, h_∞ and n_∞).

x	αx	βx	x∞
m
h
n

The leak current for Cl^- has a constant membrane conductance (g_Cl), as do those of Na⁺ and K⁺ (g_Na^l and g_K^l); E_eq was the resting potential, set here to be −70 mV.

The membrane potential (E_m) was calculated by applying Kirchoff’s law with the assumption that the cell membrane acts as a capacitor, with C_m denoting membrane capacitance: (7)

The Nernst potentials (v_ion) of an ion with valence z_ion are given by: (8) where R is the universal gas constant and T is the absolute temperature.

The pump current rate was modelled by the formula introduced by Luo and Rudy [18]: (9) where P_max represents the maximum pump current, and K_m, Na_i and K_m, K_e the half saturation concentrations for [Na]_i and [K]_e, respectively. f_NaK reflects the effect of rectification of membrane potential on the Na/K pump: (10)

Modelling cell volume

The volume of a cell (v_c) was assumed to be regulated in response to the osmotic pressure difference, which drives water into or out of the cell: (11) where L_p is the membrane permeability to water flow, and A_s is the cell membrane surface area. The cell membrane surface area was assumed here to be constant for simplicity despite changes in the shape and hence the volume of the cell. Non-permeable organic molecules inside cells also contribute to cell volume regulation, and their concentration (A_i) varies inversely with intracellular volume, assuming a fixed total amount present.

The interstitium was assumed to be a continuum, in which cells are uniformly distributed at constant number density (N_c) and collectively constitute the intracellular space, surrounded by the extracellular space. Making the assumption of homogeneous cells, the intracellular volume at tissue scale is given by V_i = v_c N_c. The volumes of intracellular and extracellular space were allowed to vary according to V_i + V_e = V₀, where V₀ is a representative elementary volume. In terms of volume fraction, the relationships read: (12) (13) where ∅_i and ∅_e are the volume fractions of intracellular and extracellular space, respectively.

Modelling ion transport in extracellular space

In addition to the local influx and efflux of ions across the cell membrane locally, ions also transport via the mode of diffusion (concentration gradient as the dominant driving force, neglecting electrodiffusion) in the extracellular space. Extracellular ion concentrations were derived from the mass balance of the ion species: (14) (15) (16) where D_ion^* are the effective diffusion coefficients of the ions in the extracellular space, incorporating the effect of varying extracellular volume fraction on diffusion coefficients in free media as: (17)

The term S_ion represents a source/sink term for ions, i.e. influx/efflux of ions across the cell membrane is reflected at the tissue scale with the assumption of a collection of homogeneous cells lumped into the intracellular space: (18) (19) (20)

In the current study, the geometry has been set to be a one-dimensional spherically symmetrical coordinate space. Therefore, conservation equations of extracellular ions can be expressed explicitly as: (21) (22) (23)

Modelling thermal cell death

It is of great importance to have proper descriptions relating tissue temperature to tissue state; this was chosen as cell viability in the current study. The hyperthermic cell death model developed by O’Neill et al. [11] was adopted here, since this takes into account both fast cell death during ablation (cell transition among three states proposed as alive (denoted by A), vulnerable (V), and dead states (D)), and slow cell death days after thermal treatment. A schematic of this model is shown in Figure 1. Using standard reaction kinetics, the system is described by two simultaneous differential equations, given the constraint that A + D + V = 1: (24) (25) (26) where k_f is the forward rate whilst k_b is the backward rate. It was noted that k_f is a function of temperature, with k_f being a scaling factor and T_k a parameter that calibrates the exponential dependence on temperature; this term also incorporates the factor (1 - A) to acknowledge the potentially greater effect of the damaged state of surrounding cells on the reaction dynamics. A continuous decrease in cell viability through the slow cell death mechanism was assumed to initiate 2 h post-treatment, when cell viability was below a critical post-heating cell viability threshold (1 - D_τ). (27) (28) where k_s is the transition rate of a secondary, slow cell death mechanism and includes a baseline scaling factor k_s, and the threshold value. The incorporation of D(1 − D) was done to give the best fit to the data while satisfying the restrictions imposed on k_s (see O’Neill et al. [11] for more details).

The final part of the model was the inclusion of the impact of thermal cell death on the other sub-models.

Effect on trans-membrane ionic kinetics

It was assumed here that the cell membrane does not rupture during ablation, thus providing a continuous barrier between ion concentrations in the intracellular and extracellular spaces. It was further assumed that active ionic transport (i_NaK) was disrupted when cells transition into the dead state, while having a negligible effect on passive ionic transport (i_Na, i_K, and i_Cl). However, passive ionic transport was also affected as it was essentially determined by the ion concentrations. In other words, upon cell death, ionic transmembrane transport was purely a passive process, dependent on chemical and electrical driving forces without biological regulation through Na-K pumps.

Despite there being heterogeneous cell types (dead versus alive) induced by thermal treatments, ‘averaged’ intracellular ion concentrations at the tissue scale can be approximated by linearly correlating the fraction of cell viability to the active ionic currents. Note that the vulnerable (V) cell state is the cell state located between the fully alive (A) and dead (D) states. In a biological sense, it is intended to represent a state in which the cell’s normal functions may be impaired but when the cell is not yet fully dead. It has been assumed here that only fully alive cells are capable of powering ion pumps and regulating osmotic volume to maintain homeostasis during and after thermal ablation. Thus, the modified form of intracellular Na⁺ mass balance at the tissue scale (Equation 18(18) ) is given as: (29)

Similarly for K⁺, the final expression is: (30)

Effect on cell volume

It was assumed that the change of cell volume undergoes two phases. In the first phase the cell volume follows osmotic regulation during ablation, and then it is deprived of such regulation due to non-viability of cells within the ablated zone, upon which it stops changing and remains at the final volume independent of further change in transmembrane ion concentrations. Thus Equation 12(12) describing change of volume fraction at the tissue scale was modified by correlating it with cell viability; this permitted a consistent description of cell volume change by osmosis in the unablated region: (31)

Numerical methods

In total, the model contains 10 independent variables, which can be written in vector form as u = ([Na]_e, [K]_e, [Cl]_e, [Na]_i, [K]_i, [Cl]_i, ∅_i, E_m, A, D)^T. The system of governing equations was correspondingly finalised by Equations 21(21) , 22(22) , 23(23) , 29(29) , 30(30) , 20(20) , 31(31) , 7(7) , 24(24) and 25(25) /27(27) . Three of the governing equations are partial differential equations dealing with ion diffusion in the extracellular space while the others are intrinsically ordinary differential equations (ODE) based. When the space-dependent terms are discretised in a standard finite difference method, the spatial derivatives are converted into a source term, which allows the system of equations to be written in the form (32) where M corresponds to the mass matrix in accordance with the temporal derivative vector of the 10 variables, whilst F_diff are forcing functions associated with extracellular diffusion. The resulting system of ODEs was solved in MATLAB using a stiff solver, with ode15s being chosen based on considerations of numerical accuracy as well as the time and state dependence of the mass matrix.

Temperature function

For simplicity, a representative tempo-spatial temperature profile during ablation [19] was assigned to initiate simulations: (33) where T₀ is the baseline temperature and ΔT is the maximum temperature rise. The hyperbolic tangent function is specified by parameters (r_c, t_c,1 and t_c,2) that set the offsets of the centre points, and parameters (r_s, t_s,1 and t_s,2) that scale the stretches of spatial and temporal distributions. The temperature profile is depicted in Figure 2 using the given baseline parameter values.

Figure 2. Temperature profile imposed on the system.

Parameter values

Parameters in the study were roughly divided into three categories: those involved in the ion-based model, those relating to ion transport, and those describing thermal cell death. For the typical Hodgkin–Huxley type of ionic currents in the present study, values of the conductance parameters were chosen to ensure that the system was at its steady state at t = 0, i.e. null currents of each ion species, according to initial conditions, and capable of displaying long-term stability even under physiological conditions, i.e. the variables return to their original steady state after minor perturbations to them. The initial concentrations of intracellular non-permeable anions were calculated to fulfil osmotic balance at the initial state. Diffusion coefficients of ions were used as those estimated in brain tissue, which also acts as a porous medium. Such a choice was capable of reflecting the timescale of diffusion transport in the tumour extracellular space. Finally, parameter values in Φ the thermal death model were readily available from the work by O’Neill et al. [11]. A complete list of parameter values and their sources is given in Table 2.

Table 2. Parameter values.

Name	Value and unit	Description	Reference
R	8314 J/kmol/K	Universal gas constant
e	1.6*10^−19 C	Elementary charge
F	96485 C/mol	Faraday’s constant
T0	310 K	Physiological body temperature
Cm	47pF	Membrane capacitance	[19]
As	7.67*10^−5 cm^2	Membrane surface area	[16]
vc^0	3.8*10^−8 cm^3	Initial volume of cell	[16]
Nc	0.2105*10^8 1/cm^3	Tumour cell density
Ai^0	80.35 mM	Initial concentration of intracellular non-permeable molecules	Calculated*
gNa^l	4.21*10^−6 mC/V/min	Leak conductance of Na^+	Calculated*
gK^l	7*10^−6 mC/V/min	Leak conductance of K^+	Calculated*
gNa^g	8.174*10^−4 mC/V/min	Maximum gated conductance of Na^+	Calculated*
gK^g	9.123*10^−4 mC/V/min	Maximum gated conductance of K^+	Calculated*
gCl	2.81*10^−4 mC/V/min	Leak conductance of Cl^-	Calculated*
Km, Nai	10 mM	[Na]i at half maximum pump current	[17]
Km, Ke	1.5 mM	[K]e at half maximum pump current	[17]
Pmax	73.97pA	Max. pump current	[17]
Lp	10^−9 m/Pa/s	Membrane permeability to water flow	[16]
	10.5 cm^2/min	Na^+ diffusion coefficient	[12]
	11.76 cm^2/min	K^+ diffusion coefficient	[12]
	12.18 cm^2/min	Cl^− diffusion coefficient	[12]
kf	3.33*10^−3 1/s	Scaling constant for forward rate	[11]
kb	7.77*10^−3 1/s	Backward rate constant	[11]
Tk	40.5 °C	Calibration temperature	[11]
ks	0.316*10^−3 1/s	Scaling factor for slow-death rate	[11]
Dτ	0.208	Cell viability threshold	[11]
L	3 cm	Length of computational domain	[18]

Initial and boundary conditions

Initial conditions for all variables were set to be those values determined under physiological temperature and are listed in Table 3. To complete the diffusion–convection–reaction equations of extracellular ions, a zero flux boundary condition was assigned at the origin, and baseline physiological concentrations were set at the outer radius which was assumed to be sufficiently far away from the origin.

Table 3. Initial values of variables.

Name	Value and unit	Description	Reference
[Na]e	120 mM	Extracellular Na^+ concentration	[20]
[Na]i	27 mM	Intracellular Na^+ concentration	[20]
[K]e	4 mM	Extracellular K^+ concentration	[20]
[K]i	130.99 mM	Intracellular K^+ concentration	[20]
[Cl]e	124 mM	Extracellular Cl^− concentration	[20]
[Cl]i	9.66 mM	Intracellular Cl^− concentration	[20]
[∅]i	0.8	Intracellular volume fraction	[21]
Em	−68 mV	Membrane potential	Calculated*
A	0.9999	Fraction of cells alive	[11]
D	0	Fractions of dead cells	[1]

Results

Figure 3 shows the simulation results of cell viability based on the fast cell death kinetics, which was sufficiently accurate to represent the ablation volume in response to the imposed temperature profile here. Irreversible cell death was induced and completed rapidly within minutes, generating an ablation region with zero cell viability, beyond which cells stayed undamaged.

Figure 3. Cell viability in response to the imposed temperature profile.

Figure 4 shows the spatial distribution of ion concentrations using the given temperature profile at varying time points. Shortly after the ablation, large changes in the intracellular/extracellular ion distributions were seen in the ablation region, as the activity of the Na-K-ATP pumps diminishes with decreasing cell viability. These changes are seen as the accumulation of sodium in the intracellular space with a corresponding depletion of potassium, with the reverse trend in the corresponding ions in extracellular space. Consequently, the cell membrane potential undergoes a sudden change, from a negative intracellular charge to a positive one, which in turn causes an influx of chloride ions into the tumour cells. As time goes on, deviations of the intracellular/extracellular ion concentrations from their baseline values were observed beyond the ablation region. The steep concentration gradient in the periphery of the cell death domain provides a powerful driving force for ion diffusion in the remaining domain, where cells were still alive and thus capable of regulating the intracellular ions concentrations with the functioning Na-K-ATP pumps.

Figure 4. Simulation results of ion concentrations: top panels Na_e and Na_i; middle panels K_e and K_i; bottom panels Cl_e and Cl_i.

Figure 5 shows the spatial trend in φ_i (the intracellular volume fraction). φ_i increasing slightly until cells were completely killed, then staying at a constant level in the ablation zone (from the origin to r ≈ 1.2 cm). Beyond the ablation region, φ_i was regulated in accordance with the tempo-spatial distributions of intracellular/extracellular ion concentrations. It was shown that φ_i rises rapidly to a peak and then decays rapidly, generating a peripheral rim surrounding the ablation zone with a significant increase in φ_i. In the far region, only a minor increase in φ_i was found, with a time frame of days, due to the long timescale of ion diffusion.

Figure 5. Simulation result of volume fraction of intracellular space.

The influence of the parameters defining the temperature function, including the maximum temperature rise (ΔT), temporal and spatial ‘stretch’ (t_s and r_s) of the hyperbolic tangent temperature function upon the ablation volume and φ_i, are shown in Figure 6. The radial width for which φ_i ≥90% ≥ 90% at t = 24 h t = 24h post-ablation (i.e. a 12.5% increase relative to the baseline level), denoted by r_90%, was arbitrarily chosen to reflect the rim of cell swelling.

Figure 6. Dependence of the ablation volume and r_90% upon variations in: (a) ΔT; (b) t_s; and (c) r_s.

Figure 6(a) shows that the width of the ablation region increases slightly with increasing ΔT, while r_90% shows only marginal variation with various values of ΔT (r_90% being approximately 10% of the size of ablation zones here). It is shown in Figure 6(b) that neither the ablation volume nor r_90% was dependent on t_s. Although varying t_s affects the rate of temperature rise from T₀ to the peak (T₀ + ΔT) by ‘squeezing’ or ‘stretching’ the temperature profile, thus delaying or advancing initiation of cell death, its effect was confined within the same heating zone from the origin up to r≈r_c, resulting in an identical ablation volume. As shown in Figure 6(c), the ablation volume increases with greater r_s, which spatially stretches the temperature profile, rendering a broader region with temperature rise to certain extent. However, the rim of cell swelling seems to be of approximately constant width irrespective of r_s except at r_s = 0.1.

Presented in Figure 7 is the impact of ion diffusion coefficients on the rim width, which displays a non-trivial dependence on diffusion coefficient values when they are varied in the range from one tenth to ten times baseline values. When diffusivity was reduced by a factor of 10, a wider swelling rim was found for K⁺, while the rim was narrower for Na⁺ and Cl⁻, this effect being more significant for Na⁺. On the other hand, when diffusivity was increased 10-fold, the rim width was unchanged for changes in K⁺ and Cl⁻, while a slight increase was found for Na⁺. Overall, therefore, an increase in diffusivity of Na⁺ can lead to an increase in the rim width, while the rim width is essentially independent of the diffusivity of K⁺ and Cl⁻ within a certain range.

Figure 7. Dependence of r_90% upon ion diffusion coefficients.

Figure 8 shows the spatial dependence of the intracellular cell volume fraction on various initial intracellular volume fractions. Given that the volume of a single cell is fixed, the effect of initial intracellular volume fraction is analogous to that of cell packing density. For the lower value of ∅_i0, despite a significant increase in ∅_i at the peripheral region, its intensity was somewhat weaker than that for the greater value of ∅_i at 24 h. As time progressed, the intracellular volume fraction dropped quickly for φ_{i, 0} = 0.6 = 0.6. This reflects the effect of cell packing density on ion diffusivity, which is inversely dependent on extracellular volume fraction (i.e. the increased intracellular volume fraction tends to block ion transport space).

Figure 8. Dependence of intracellular volume fraction on its initial level: (a) φ_{i, 0} = 0.8; and (b) φ_{i, 0} = 0.6.

Discussion and conclusions

In this study a novel mathematical model was developed to describe the enhanced rim on the periphery of the ablation region, which is ascribed to cellular oedema as a physiologically osmotic response. It captures the information flow, which is the redistribution of ions, initiated with ion exchange between intracellular and extracellular space locally perturbed by thermal injury, followed by global ion movement between the ablated region and its surrounding normal tissue, with cell swelling as a final consequence. As an initial effort, the mathematical model incorporates minimally essential elements in a spatially distributed system, namely components of (1) ion concentration dynamics, (2) cell volume regulation, (3) ion transport, and (4) thermal cell death. It should be noted this represents a proof-of-concept example to illustrate a potentially indirect response to thermal injury, thus supplementing evaluation of the ablation volume. Indeed, an analogous but more evident example is the immunogenic response post ablation, which has been indicated by a range of pre-clinical and clinical evidence [20–23] and which is believed to be induced by release of inflammatory cytokines from the ablated tumour cells and subsequent infiltration, and thus is observed at the margin of the ablation region as well as in distant, untreated regions post-ablation (within a time frame ranging from hours to days).

Thermal cell death generates an ablation zone, distinguishing the tissue state from ‘alive’ to ‘dead’, which indicates that ion dynamics should be dependent on the tissue state. Continuous descriptions of ion dynamics were adopted in distinct regions from the ablation zone to the undamaged regions by assuming that the effect of ablation on ion dynamics is modelled by correlating cell viability with activity of the Na-K pumps, as it has been suggested that Na-K pumps stop functioning in the temperature range for ablation [24]. Another likely biological response that could be proposed is cell membrane rupture; however, from a mathematical modelling point of view, this is computationally highly demanding due to its discontinuous nature. Were it to happen, extracellular ion concentrations would need to be corrected to account for all ions suddenly appearing in the extracellular volume; although this would be in a volume fraction averaged form, this would yield a sudden decrease in Na_e and Cl_e and a sudden increase in K_e. The same changing trend in ion concentrations was qualitatively predicted by the continuous description used here.

The intracellular volume fraction was the target output of the model, and this was regulated in response to osmosis. Regarding the choices of ion species involved and their transport modes through ion channels, the model was formulated to have general and well-established descriptions of Na⁺, K⁺ and Cl⁻, to reflect a qualitative trend in intracellular volume fraction when subject to a disturbance of ion concentrations of all species upon thermal ablation. In addition, it was assumed that intracellular volume regulation was not limited due to neglecting the mechanical force exerted by the tumour extracellular matrix, which was composed of mainly collagen fibre bundles and proteoglycan filaments, thus providing resistance to large changes in the extracellular volume. Such aspects can be incorporated preliminarily by including a linear reaction term (representing the restoring force) [16], which may weaken the extent of cell swelling; however, this is not likely to alter the qualitative trend. Furthermore, the impact of cell mechanics on cell volume regulation has not been taken into account in the current study. It is reported that changes in membrane tension tend to induce mechano-sensitive channels (passive membrane proteins) to open, through which ion transport is driven by the concentration gradient [25]. Such aspects of the behaviour could be naturally included in the current model setting and will be examined in the future.

Ion transport in extracellular matrix is the only means to provide a spatial link to connect the locally dependent modelling components of ion dynamics and cell volume regulation in a spatially distributed system. In the current model, diffusion was assumed as the only mode of ion transport without accounting for the other transport modes via convection and migration. Despite this, it is still complex to understand the highly non-linear dependence of the rim width upon varying ion diffusion coefficients, as the change in the diffusion coefficient of individual ions can lead to changes in all ion concentrations simultaneously in a highly non-linear manner, which in turn work in concert to regulate the cell volume. However, there are several areas for potential refinement of ion transport. Firstly, though the interstitial fluid velocity has been found to be negligible in the tumour core region due to uniformly elevated interstitial fluid pressure in this region it cannot be so easily neglected in the tumour boundary [26], where convection may overwhelm diffusion, driving outward ion movement from the ablation region to the periphery, thus affecting the spatio-temporal distribution of ions, and consequently the timescale and spatial extent of cell swelling.

In addition, when considering the tumour microenvironment as containing fixed interstitial proteins with negative charge (e.g. glycosaminoglycans) as well as potentially ‘ablation-enhanced’ diffusion, local built-in electric fields (self-fields) can be induced, which indicates that ion migration along voltage gradients should be included in further study. It would be of great importance to test our hypothesis and to further assess the timescale and extent of cell swelling, thus reaching a qualitatively solid conclusion of the occurrence of cell swelling and revealing any correlation between it and the ablation volume.

A mathematical model has been developed with essential descriptions of the minimally necessary components to capture cellular swelling at the margin of the ablation region after thermal treatments, providing a possible cause of the bright rim observed from imaging and a valuable link between the interpretation of the underlying physiological process and post-treatment imaging. In addition to the primary response to the thermal injury the model demonstrates a potential interaction between the ablated cells and the non-ablated cells. It is worth noting that one potential interaction that the upstream hyperthermic zone has is cell detachment in the downstream non-hyperthermic zone; this has been demonstrated in an in vitro microchannel [27], which was illuminating within our context for experimental verification of cell swelling within a controllable experimental platform. Alongside experimental progress and advancing imaging techniques, the model could also be of great use in two ways: on one hand the model serves as a simulation platform and is capable of incorporating more detailed descriptions (i.e. ion transport), on the other hand it could be integrated into the already validated mathematical models for thermal therapies (e.g. a multi-scale model for radiofrequency ablation in liver tumours by Payne et al. [28]), which provides both the outputs of the bio-heat model and of the cell death model with a known degree of confidence, thus allowing for the testing of predictions of the cell swelling rim in a more realistic scenario. In a clinical setting, the location of this bright rim could be of considerable potential use for determination of the extent of the ablation region.

Disclosure statement

The authors alone are responsible for the content and writing of the paper.

Funding

The research has been funded by the European Commission under [Grant Agreement No. 600641], FP7 Project Go-Smart.