1. Introduction

Organ transplantation remains one of the most effective therapies for end-stage heart, kidney, and liver failure, offering patients the possibility of functional recovery and prolonged survival^[1-3]. Since the first successful kidney transplantation in 1954, the field has expanded rapidly, with transplantation now performed at a large scale worldwide^[4]. However, ischemic injury begins to accumulate immediately after organ procurement, resulting in a narrow time window for maintaining graft viability^[5]. Donor hearts preserved by conventional static cold storage (SCS) are generally optimal when ischemic time is kept below approximately 4-6 h, beyond which the risk of graft dysfunction increases markedly, contributing to organ discard and exacerbating the organ shortage^[6,7].

In current clinical practice, SCS remains the prevailing preservation approach and has long been regarded as the clinical standard^[8]. Its preservation time, however, is limited for metabolically active organs, particularly the heart and liver. Machine perfusion (MP) has therefore attracted increasing interest because it can extend preservation duration while mitigating hypoxic injury by maintaining controlled perfusion and enabling oxygen delivery. Evidence from clinical and experimental studies indicates that MP can improve preservation outcomes compared with SCS in selected settings. Most established preservation strategies operate at or above 0 °C (typically 0-4 °C), whereas subzero organ preservation remains less explored^[9]. In principle, further temperature reduction can prolong viable storage and improve preservation outcomes^[10]. Conventional SCS at subzero temperatures leads to ice formation and severe tissue injury, rendering grafts unsuitable for transplantation. Supercooled organ preservation addresses this limitation by delivering cryoprotective agents (CPAs) via MP and by precisely controlling cooling and rewarming, enabling the organ and perfusate to remain ice-free below 0 °C. Therefore, supercooling is promising for extending preservation time because metabolic activity decreases substantially with temperature; a widely used approximation is that metabolic rate decreases by 50% for each 10 °C reduction^[11].

Supercooling is intrinsically metastable. Maintaining an ice-free state requires stringent temperature stability because thermal perturbations can trigger intracellular ice formation and cause lethal cellular injury^[12]. As tissue scale increases, limited mass transport, pronounced heterogeneity, and slow heat transfer become increasingly prominent, making the suppression of heterogeneous ice nucleation within large organs progressively more difficult^[13]. At the same time, constrained by current temperature-sensing techniques and ice-crystal observation methods, it remains difficult to monitor intra-organ temperature fields in real time and to proactively assess freezing risk during supercooled preservation^[14]. In addition, organ-specific anatomical features and vascular architectures lead to spatially varying perfusion, which can create localized high-risk zones for freezing^[15].

To address these challenges, researchers have developed bioheat models and numerical simulation frameworks to predict intra-organ temperature distributions and perfusion-rate fields during preservation^[16-18]. In silico prediction of thermal behavior can identify regions at risk of freezing, while quantitative evaluation of perfusion dynamics can guide protocol design, reduce experimental burden, and improve preservation outcomes. Berendsen et al.^[19] combined supercooling with MP in rat livers and achieved long-term transplant survival after 4 days of preservation. Ozgur et al.^[20] optimized subzero preservation workflows and reported that low-flow reperfusion during recovery can benefit graft outcomes. Quader et al.^[21] investigated ex vivo perfusion rates matched to cardiac metabolic demands across temperatures and identified flow conditions required to sustain aerobic metabolism. These studies highlight the potential of supercooled MP and indicate that resolving perfusion-rate distributions is important for achieving favorable preservation performance. Because direct measurements of intra-parenchymal temperature and perfusate hydrodynamics remain limited, computational bioheat modeling is indispensable for mechanistic analysis. Finite-element simulations have been used to resolve evolving tissue temperature fields during perfusion and to enable parametric sweeps for evaluating thermal outcomes and optimizing protocols. Deng et al.^[22] simplified renal vasculature into straightened vessels based on Murray’s law and assessed how flow and pressure vary with vessel number and branching patterns. Emerson et al.^[18] introduced a space-colonization algorithm and the “super-lobule” concept to model hepatic vascular flow and identify poorly perfused regions. However, these studies still rely heavily on algorithmically generated idealized geometries; the irregular morphology and complex vascular trees of real organs can produce substantial deviations in temperature and flow fields relative to ideal models. Consequently, anatomically based perfusion–heat-transfer models have gained momentum. Ehrlich et al.^[23] performed thermal analysis of vitrification-based preservation using rabbit-kidney data to inform human-kidney modeling, although experimental validation remained limited. Moreover, while supercooled perfusion has been demonstrated in small-animal organs, translation to large-animal and human-scale organs remains to be established.

In this study, we leveraged a self-developed supercooled MP platform to investigate thermal effects during multi-organ supercooled perfusion. We developed a multi-organ model for the liver, heart, and kidney and predicted intraparenchymal temperature distribution under the condition of realistic supercooled machine perfusion operation. Model fidelity is assessed by comparing simulated temperatures with multi-site thermocouple measurements. This work systematically integrates organ inlet perfusion parameters, convective heat transfer at the organ surface, perfusate thermophysical properties, and a multi-organ physical model to elucidate the influence of these factors on the heat transfer processes in supercooling MP. The study provides theoretical guidance and a reference for thermodynamical optimization in supercooling MP protocols.

2. Methods

2.1 Design and assembly of the multi-organ supercooling perfusion system

We developed a multi-organ supercooled machine-perfusion system comprising a perfusate circulation loop and a refrigeration loop, in which forced convection heat transfer between the two loops ensures uniform temperature distribution and precise thermal regulation to enable ice-free preservation of donor organs below 0 °C (Figure 1). Ethanol was selected as the circulating coolant in the refrigeration loop to ensure reliable operation below 0 °C and low fluid viscosity at negative temperatures. The refrigeration unit consists of a plate heat exchanger, compressor, condenser, and centrifugal pump, with cooling capacity automatically modulated by a programmable logic controller (PLC) according to predefined perfusion temperature curves.

Display Full Size

Figure 1. Schematic of the multi-organ supercooled perfusion system and key research components. Overall layout includes supercooling organ preservation and parametric study of the heat-transfer model. Created in BioRender.com. Q: inlet flow rate; P_in: inlet pressure; h: effective external convective heat-transfer coefficient; c_p: perfusate specific heat capacity; k: perfusate thermal conductivity; ΔT = T_measure - T_simulation.

Briefly, each preservation experiment followed a standardized workflow. Freshly procured organs were first immersed in precooled University of Wisconsin (UW) solution for initial thermal equilibration, then cannulated and connected to the organ-specific perfusion circuit. After the circuit was primed and the target perfusion pressure or flow rate was established, cryoprotective agents were introduced by gradient loading, and the organs were cooled in a controlled manner to the target supercooled temperature range (typically -4 to -6 °C). Throughout perfusion, perfusate and intraparenchymal temperatures, as well as inlet and outlet pressure and flow rate, were continuously recorded until completion of the preservation protocol.

Cannulation strategies differ by organ: (1) For the heart, perfusate is delivered through the aorta and exits via the pulmonary artery; (2) For the liver, both dual and single perfusion configurations are available. In MP of liver grafts, single-vessel perfusion via the portal vein simplifies cannulation and reduces procedural risk. However, because the hepatic artery provides a substantial portion of the liver’s oxygen supply, single-vessel approaches may deliver less optimal oxygenation than dual perfusion strategies^[24]. Experiments suggest that dual-vessel perfusion via both the portal vein and the hepatic artery more closely approximates physiological hepatic blood supply^[25]. Therefore, because supercooled perfusion markedly reduces metabolic activity and oxygen demand at hypothermic temperatures, a single-vessel approach can be sufficient and advantageous; (3) For the kidney, perfusate enters through the renal artery and exits via the renal vein. The organ chamber is submerged within a coolant bath connected to the refrigeration loop and accommodates organs of varying sizes. To mitigate freezing risk during preservation, CPAs were introduced via gradient loading under MP and combined with controlled cooling-rewarming ramps to reduce ice-nucleation propensity. Rapid temperature fluctuations and localized cold spots were avoided to minimize the probability of spontaneous nucleation, and the perfusate formulation was optimized by incorporating polyethylene glycol, trehalose, and glycerol to reduce intra- and extracellular ice formation risk^[26]. Given organ-specific differences in hemodynamic tolerance, perfusion parameters were tailored to each organ to avoid pressure-induced injury. The peristaltic pump operated in closed-loop control to drive perfusate flow, while feedback from pressure sensors and flowmeters was used to regulate organ-specific perfusion pressure and flow rate. The control module implemented a controlled cooling rate to decrease perfusate temperature to the target supercooled range, typically -4 to -6 °C, thereby limiting rapid thermal excursions that could initiate ice nucleation. Throughout perfusion, flow rate and perfusate temperature were adjusted in response to temperature-probe feedback to maintain a uniform core temperature within the organ.

To monitor intra-organ temperature and detect potential crystallization, two thermocouples were used in each experiment^[27,28]. One T-type thermocouple probe was positioned within the organ chamber near the organ inlet to measure perfusate temperature, while the second thermocouple was inserted approximately 20 mm below the organ surface into the organ parenchyma to record local tissue temperature in real time. Temperature signals were continuously acquired using a high-precision data acquisition system (DAQ970A, Keysight Technologies, USA) and recorded with the accompanying software. Measurements of organ parenchyma were obtained from three independent experiments for each organ type (n = 3). The placement of thermocouples inserted into the organ parenchyma was determined based on previous studies to ensure representative and reliable temperature measurements. In the liver, the probe was positioned in the right lobes to capture regional temperature variation^[29]. In the heart, the thermocouple was placed within the myocardium of the left ventricular walls^[30]. In the kidney, the probe was inserted into the cortical regions^[31]. The specific placement site of the thermocouple in each organ is shown in Figure S2. These sensors provided real-time measurements of both perfusate and local tissue temperatures during perfusion. During perfusion, the system also monitored inlet and outlet pressure and flow rate to document relevant perfusion parameters. No ice formation was observed in either the organ tissue or the perfusate throughout the entire procedure, confirming successful maintenance of the supercooled state. The protocol and instrumentation were repeated for each organ type, with cannulation configurations and flow settings adjusted to accommodate organ-specific anatomical and physiological differences.

Overall, the integrated thermal-fluidic control and real-time perfusate and intraparenchymal temperature monitoring enabled reproducible, ice-free supercooled MP with uniform organ core temperatures across heart, liver, and kidney preservation runs.

2.2 Geometric model of multiple organs

A finite-element model was implemented in COMSOL Multiphysics 6.3 to supplement the experimental measurements and to elucidate intra-organ heat-transfer mechanisms during supercooled MP. The model resolves convective heat transfer associated with perfusate transport through the organ vasculature and conductive heat transfer within the parenchyma. The three-dimensional heart, kidney, and liver geometries were obtained from an anatomically based dataset released through the BodyParts3D database, which was constructed from whole-body magnetic resonance imaging acquired at 2 mm slice intervals from an adult male volunteer and subsequently segmented to represent organ morphology^[32]. Prior to import into COMSOL, the surface meshes were edited and retopologized in Blender. Because the available anatomical models include only major vessels and do not provide small arteries and capillaries below the imaging resolution, the parenchyma was represented as a fluid-saturated porous continuum, whereas the large vessels were modeled as free-flow lumens.

From a structural perspective, organ parenchyma comprises muscle tissue, connective tissue, and adipose components that form an interconnected three-dimensional pore network with spatially heterogeneous pore sizes and permeability^[33]. Biological tissue transport has therefore been widely formulated within porous-media theory to capture coupled flow and heat transfer at the organ scale^[34,35]. On this basis, the present model treats intra-organ heat and mass transport within a porous-medium framework, while retaining explicit free-flow regions for the major vasculature. To render the supercooled machine-perfusion process computationally tractable, the following assumptions were adopted:

(1) According to the previous bioheat and thermal-tissue modeling study^[36], material properties were treated as homogeneous and isotropic within each domain.

(2) Due to the addition of low freezing point CPAs, continuous thermocouple measurements showed no signatures of phase change. As a result, macroscopic freezing did not occur under the present experimental conditions, supporting the assumption that latent heat effects were neglected.

(3) Metabolic heat generation was neglected due to strongly suppressed metabolic activity at hypothermic and subzero temperatures.

(4) The perfusate in the organ chamber was assumed to be well mixed, and the organ surface was modeled using a convective heat-transfer boundary condition.

(5) Organ deformation was neglected, and the perfusate was treated as incompressible.

(6) Due to the macroscopic thermal diffusion time T_sys = L²/α being significantly larger than the local solid-fluid thermal exchange time T_ex = ρc_pV/h_intA, local thermal equilibrium (LTE) was assumed in the parenchyma, meaning that the tissue solid phase and the perfusing fluid are at the same local temperature.

2.3 Bioheat transfer model of organ as porous media

We solved a conjugate heat-transfer problem spanning the vascular lumens and the organ porous parenchyma. In COMSOL Multiphysics 6.3, the heat transfer in fluids interface was used to model thermal transport in the flowing perfusate within the vascular lumen and in the organ porous parenchyma, while the heat transfer in porous media was used to model heat transfer within the organ parenchyma. Within the porous-medium formulation, COMSOL provides two commonly used representations of the energy balance: a LTE model that uses a single temperature field shared by the fluid and porous matrix, and a local thermal nonequilibrium model that resolves separate temperatures for the solid matrix and the fluid phase^[37]. Under the low-flow operating conditions of supercooled perfusion, the parenchymal tissue temperature is expected to remain close to the perfusate temperature, and the LTE model was therefore adopted for the present simulations. All physical interfaces are coupled using the Nonisothermal Flow multiphysics coupling, ensuring that temperature-dependent material properties are consistently shared between the momentum and energy equations. Therefore, the heat-transfer equation for the unconstrained fluid is applied within the vasculature, and it is defined as follows:

(1)${\rho }_{\mathrm{ff}}{c}_{p,\mathrm{ff}}(\frac{\partial T}{\partial t}+u_{\mathrm{ff}}\cdot \nabla T)=\nabla \left(k_{\mathrm{ff}}\nabla T\right)$

where ρ_ff is the density of the free fluid; c_p,ff is the specific heat of the free fluid; k_ff is the thermal conductivity of the free fluid; u_ff is the flow velocity in the free flow.

For the organ parenchyma, heat transfer is modeled using the porous-media heat-transfer equation:

(2)${\left(\rho c_p\right)}_{\mathrm{eff}}\frac{\partial T}{\partial t}+{\rho }_{\mathrm{pf}}{c}_{p,\mathrm{pf}}{u}_{\mathrm{pf}}\cdot \nabla T+\nabla \cdot q=Q_h$

where (ρc_p)_eff is the effective heat capacity of porous media, which is a combination of the fluid and solid components: (ρc_p)_eff = (1 - ε)ρ_psc_p,ps + ερ_pfc_p,pf; c_p,pf is the specific heat capacity of the flushing fluid; c_p,ps is the specific heat capacity of the solid component; ε is the porosity of the organ parenchyma; u_pf is the flow velocity in the porous media; q is the heat flux vector: q = -k_eff∇T; Q_h is the general volumetric heat source.

2.4 Flow model of free flow and porous media

In this section, we describe the development of a flow model for organ supercooled MP, which can be used to assess appropriate perfusion flow rates for multi-organ supercooling preservation. Perfusate flow within the organ was modeled using two distinct regimes: free flow within vascular lumens and porous-media flow in the parenchyma. Given that the perfusate velocity in vascular lumens is relatively low, free flow was modeled using a laminar flow physics model governed by the Stokes equations. In contrast, flow within the organ parenchyma is treated as flow through a porous medium because the tissue can be approximated as a saturated, single-phase porous domain. For this porous-media region, we employ Darcy’s law to describe perfusate transport, assuming low Reynolds numbers and negligible inertial effects, consistent with laminar, permeability-dominated flow. Accordingly, turbulence effects are neglected in the porous regions due to the small local velocities.

Free flow in vascular lumen momentum and continuity equation:

(3)$\{\begin{array}{c}{\rho }_{\mathrm{ff}}\frac{\partial u_{\mathrm{ff}}}{\partial t}+{\rho }_{\mathrm{ff}}(u_{\mathrm{ff}}\cdot \nabla )u_{\mathrm{ff}}=\nabla \cdot (-p_{\mathrm{ff}}I+\tau )\\ \nabla \cdot u_{\mathrm{ff}}=0\end{array}$

where p_ff is the pressure of the free flow; I is the identity tensor; τ is the viscous stress tensor: τ = μ(∇μ_ff + (∇μ_ff)^T).

Porous media (organ parenchyma) Darcy law and mass equation:

(4)$\{\begin{array}{c}\frac{\partial }{\partial t}\left(\epsilon {\rho }_{\mathrm{pf}}\right)+\nabla \cdot \left({\rho }_{\mathrm{pf}}{u}_{\mathrm{pf}}\right)=Q_m\\ u_{\mathrm{pf}}=-\frac{K}{\mu }\nabla p_{\mathrm{pf}}\end{array}$

where Q_m is the mass source term; K is the permeability of the porous media; p_pf is the pressure of the porous flow.

2.5 Boundary conditions and material properties

The flushing fluid within the organ is treated as a mixed region, so there is no need to mesh the water tank in the flushing fluid. At the same time, since the organ is immersed in the flushing fluid during the mechanical perfusion process, the boundary conditions at the surface are convective boundary conditions. Therefore, the boundary condition at the organ surface is as follows:

(5)$-k_s\nabla T\cdot n=h(T_{\text{surf}}-T_{\text{bath}}(t\left)\right)$

where the convective heat transfer coefficient h is determined according to previous studies^[38,39], as shown in Table 1; k_s is the thermal conductivity of the organ tissue; T_surf is the surface temperature of the organ; T_bath is the bath temperature.

Kidneys, livers, and hearts from healthy Bama miniature pigs were used in ex vivo subzero MP experiments. A total of 9 pigs, with a body weight range of 30~40 kg and aged 15 months, were included in this study. For each organ type, three independent organs were used for model validation experiments (n = 3 per organ type). Each organ was connected to a custom perfusion circuit in which chilled perfusate was continuously circulated. The perfusate consisted of UW preservation solution, a widely used cryoprotective medium in organ preservation research. Parameter values and organ properties are summarized in Table 1.

Fresh organs were harvested from Bama pigs (Figure S1). The organs were immediately placed in UW solution and held statically at 5-7 °C for 3 h to stabilize their initial temperature within the 5-7 °C range prior to perfusion. Concurrently, perfusion parameters for each organ were set based on values reported in the literature, and the perfusate was precooled to match the organ’s initial temperature (Table 2). For the heart, the inlet condition was specified by mean perfusion velocity, consistent with coronary-flow-governed perfusion in ex situ heart models^[41]. For the liver and kidney, the inlet condition was specified by perfusion pressure to reflect common pressure-target protocols used to mitigate perfusion-related injury^[24,42] (Table 2).

2.6 Mesh and time step independence analysis

Spatial discretization plays a decisive role in the accuracy and fidelity of numerical predictions. In practical applications, however, blindly increasing the mesh resolution can make the computational cost prohibitive, while the improvement in the quantities of interest becomes marginal once the solution approaches a sufficiently refined regime^[43]. To achieve an effective balance between accuracy and cost, higher spatial resolution is applied in regions expected to exhibit steep gradients or where the target outputs are evaluated, while a coarser discretization is maintained elsewhere^[44]. Consequently, mesh independence and time step independence studies are performed to confirm that the reported results are not artifacts of the chosen spatial and temporal discretization parameters.

Using the mesh generation tools and the time dependent solver time stepping controls in COMSOL Multiphysics, four mesh resolutions and four-time step sizes were evaluated for each organ under an identical cooling protocol. Temperature curves at selected monitoring locations were then compared to assess convergence. All mesh configurations produced highly consistent temperature versus time curves, with deviations below 1%, indicating satisfactory convergence (Figure 2a,d,g). Considering the requirements for both computational cost efficiency and solution accuracy, the mesh of 243,965, 204,615, and 387,400 polyhedral cells was selected for subsequent simulations of liver, heart and kidney respectively. Figure 2b,e,h shows that a time step of 10 s is sufficient to achieve the required accuracy. Therefore, a uniform time increment of 10 s was used in all subsequent simulations to accommodate both solution fidelity and computational cost.

Display Full Size

Figure 2. Mesh convergence and time step convergence assessment for transient organ cooling simulations. (a-c) Liver; (d-f) heart; (g-i) kidney. (a, d, g) Temperature–time histories at selected monitoring locations computed using four mesh resolutions (element number n as labeled); (b, e, h) Corresponding temperature-time curves computed using four solver time-step sizes (dt as labeled); insets show magnified views of the highlighted intervals; (c, f, i) Representative finite-element surface meshes of the liver, heart, and kidney models used in the convergence assessment.

2.7 Model validation

Experimental data are critical for the development and validation of computational models of heat transfer and flow. To assess the accuracy of the present model, simulated temperatures were compared against experimentally measured thermocouple data. Firstly, calibration was specifically performed on the cooling curve of the perfusate. Although the cooling curve of the perfusate was set via the PLC, deviations from the set curve were observed at the organ inlet due to factors such as heat loss through the tubing. To minimize this systematic error, a thermocouple was placed at the organ inlet to directly measure the perfusate temperature, thereby adjusting the inlet perfusate temperature curve in COMSOL. The boundary condition of the inlet perfusate temperature was adjusted using experimental measurements to ensure that the temperature trends of the perfusate closely matched the experimental data. After this calibration process, the model’s predicted temperature curves at the organ inlet were rigorously compared with independently measured perfusion data to validate the model’s ability to capture the thermal response of each organ. The simulated inlet perfusate temperature closely matched the experimentally measured outlet temperature curves, as shown in Figure S3. This approach follows standard practices in bioheat modeling, in which model predictions are evaluated against experimental measurements to establish credibility and confidence in the simulation results. To represent organ morphology in the numerical model, three-dimensional organ geometries were derived from the BodyParts3D repository and imported into the simulation workflow. These geometries were then discretized using unstructured finite-element meshes (Figure 2c,f,i). The simulation results were subsequently compared with thermocouple measurements recorded throughout the entire duration of each perfusion experiment (n = 3 for each organ). For all three organ models, the predicted temperature curves closely tracked the experimental curves, although minor differences were observed at certain stages of the cooling trajectories. Nevertheless, the overall trend between simulation and measurement was highly consistent. Observed deviations were documented, and their implications are discussed in Section 3. Importantly, experimental validation confirmed that the model reliably predicts the thermal response of organs under subzero perfusion conditions, thereby enhancing confidence in its use for further parametric studies and optimization of perfusion protocols. Similar validation approaches, in which numerical predictions are benchmarked against measured temperature data to establish model credibility, are commonly employed in bioheat transfer modeling studies.

3. Results and Discussion

3.1 Analysis of thermal distribution in multi-organs

To validate the thermal model, simulated organ temperature evolution curves were compared with experimental measurements from multiple thermocouples integrated into the multi-organ perfusion apparatus. Figure S4a,b,c shows the measured temperature trajectories for the liver, heart, and kidney, exhibiting consistent temporal trends across replicates. Figure 3 shows the organ temperature curves at selected monitoring locations for the liver, heart, and kidney, comparing COMSOL Multiphysics simulation results with corresponding thermocouple measurements at the same anatomical locations. In addition to bulk temperature trajectories, we also examined the spatial temperature distributions within each organ by extracting simulated temperature fields at corresponding anatomical cross-sections.

Display Full Size

Figure 3. Experimental validation of transient cooling and intra-organ temperature non-uniformity during supercooled machine perfusion of porcine organs. (a) Liver, (c) Heart and (e) Kidney parenchymal temperature time curves comparing measurements and simulations (n = 3 independent experiments, 1 organ in each group); (b, d, f) Spatial distributions of the temperature difference between intra-organ tissue and perfusate, ΔT’ = T_tissue - T_perfusate, for the (b) liver, (d) heart, and (f) kidney at the representative time points labeled in each panel. Color legends indicate ΔT’ (K).

Quantitatively, the model predictions closely tracked the experimental curves throughout the perfusion period, with only minor deviations. In most instances, the simulated and measured temperatures of the liver and kidney differed by less than 1 °C at any given time point, a discrepancy comparable to the typical accuracy of the thermocouple sensors (Figure 3a,c,e). Across all measured temperature trajectories, the model exhibited mean absolute error (MAE) ± SD of 0.24 ± 0.14 °C (n = 3) for the liver, 2.63 ± 1.02 °C (n = 3) for the heart, and 0.40 ± 0.24 °C (n = 3) for the kidney (Figure S4d,e,f). The maximum temperature differences were 0.90 °C, 5.24 °C, and 1.81 °C, respectively. This level of agreement indicates a high degree of fidelity in the model’s thermal predictions under the tested conditions. However, small differences about 5 °C between simulated and measured temperature were observed in the heart model (Figure 3c), which can be attributed to known experimental and modeling factors. From an experimental standpoint, the discrepancy between measured and simulated temperature profiles in the heart may be exacerbated by enhanced heat loss at lower perfusion temperatures. As the perfusate temperature decreases, the gradient between the cardiac tissue and the ambient environment increases, which can lead to greater environmental heat dissipation and larger deviations in the observed temperature curves. From a modeling perspective, our current framework assumes homogenized physical properties for the heart. However, unlike the liver and kidney, myocardial tissue exhibits pronounced structural and functional anisotropy^[45-47]. This simplification likely contributes to the larger errors observed in the cardiac temperature predictions. In addition, the perfusate flow paths in the heart are dominated by the coronary arterial network, which was not incorporated in the present model. The absence of explicit coronary circulation representation limits the model’s ability to capture realistic fluid flow and heat transfer dynamics, further contributing to the observed discrepancies. These findings underscore the critical role of accurate anatomical and flow pathway representation in thermal simulations of supercooled mechanical perfusion, particularly for organs with complex vasculature such as the heart. Figure 3b,d,f show the simulated spatial temperature distributions of the liver, heart and kidney at the representative time points.

Taken together, these factors explain the modest deviations, which are insufficient to substantially affect the overall comparison of cooling trends. Despite these minor differences, there is strong correspondence between the model and experimental data. The model successfully reproduced the key features of organ cooling under subzero MP, including the rapid initial temperature decline, the slowing of temperature change as the organ approaches the target perfusate temperature, and the eventual stabilization of organ temperature. Furthermore, the consistency between simulated and measured spatial temperature distributions supports the model’s capability for accurately representing intra-organ thermal distribution. For the majority of the perfusion duration, the experimental curves lie closely around the simulated curves, indicating that the computational model accurately represents the dominant heat-transfer mechanisms, convective cooling via perfusate delivery and conductive transport within the tissue. This agreement between simulation and measurement indicates that the thermal model for porous-media is validated against experimental data and demonstrates sufficient fidelity to simulate thermal behavior in the multi-organ system under the tested conditions. However, the current framework is expected to provide reliable predictions primarily for organs consistent with these assumptions, whereas its application to highly anisotropic organs without detailed vascular flow modeling should be interpreted with caution. This also suggests that vascular structure plays an important role in determining the accuracy of thermal predictions.

Overall, the simulation captured the cooling trends exhibited by all three types of organs. In each case, organ temperature declined rapidly after the onset of perfusion and then gradually approached the perfusate temperature, reaching a stable state; this characteristic behavior was reproduced in both the simulated and measured data. The simulated temperature maps further revealed coherent spatial gradients, with cooler regions propagating from major perfusion conduits into the surrounding parenchyma, consistent with expected heat-transfer mechanisms. However, larger deviations were observed for the heart compared to the liver and kidney, so the model should be used with caution for organs with more complex vascular architectures. Therefore, the model is considered suitable for predicting organ temperature evolution and spatial temperature distribution in subzero perfusion devices and provides a validated, under the tested conditions, tool for exploring and optimizing cooling protocols in multi-organ preservation.

3.2 Parameterization of perfusion inlet conditions

We performed a parameterized investigation of inlet perfusion conditions to quantify the impact of inlet settings on cooling performance in the multi-organ system. Specifically, a series of independent transient simulations was conducted in which the inlet perfusion pressure and inlet flow rate were varied within experimentally relevant ranges. For each simulation case, temperature time curves were recorded at representative monitoring locations within each organ, and the corresponding transient intra-organ temperature fields were extracted from the computed temperature distributions. The objective of this parametric study was to identify inlet settings that achieve more rapid cooling and improved temperature uniformity across organs under a recirculating perfusion configuration.

Figure 4 shows that the cooling rates of the liver, heart, and kidney in the absence of perfusion fluid are substantially lower than those observed with perfusion fluid. The simulated final-time temperature cross-sections for the liver, heart, and kidney (Figure 4b,d,f) further illustrate that higher perfusion rates enhance temperature uniformity within each organ. We introduce the maximum temperature difference ΔT_max to quantify the effect of the input perfusion condition on temperature uniformity (Figure S5a,b,c). These results demonstrate that perfusion fluid flow significantly accelerates organ cooling by improving convective heat transport. Organs with higher perfusion rates at the final time point exhibited improved temperature uniformity, although the magnitude of this benefit varied considerably among organs.

Display Full Size

Figure 4. Organ parenchymal temperature response curves under varied perfusion inlet conditions during supercooling machine perfusion. (a) Liver temperature curves over time at varied inlet pressures (0, 10, 30, 50, 70 mmHg); (c) Heart temperature curves over time at varied inlet velocity (0, 0.01, 0.05, 0.1, 1, 10 mm/s); (e) Kidney temperature curves over time at varied inlet pressures (0, 50, 100, 150, 200 mmHg); (b, d, f) Final time-point organ parenchymal temperature slices showing the spatial temperature distribution in the (b) liver, (d) heart and (f) kidney, respectively, under different inlet perfusion parameters. Color legends indicate ΔT (Κ).

The parenchymal temperature of the liver shows minimal dependence on the perfusate flow rate compared with smaller organs (Figure 4a). This observation suggests that larger biological tissues are less sensitive to perfusion-enhanced heat transfer. The differential temperature response among organs may be attributed to variations in vascular network density and tissue heat capacity, with the liver’s larger anatomical size and complex perfusion pathways limiting effective heat transfer. As shown in Figure 4a,c,e, beyond approximately 30 mmHg for liver perfusion pressure, 1 mm/s for cardiac perfusion velocity, or 50 mmHg for renal perfusion pressure, further increases in perfusion produced only negligible changes in temperature uniformity. This suggests that under high perfusion inlet conditions (> 30 mmHg for liver, > 1 mm/s for heart, 50 mmHg for kidney) heat transfer is increasingly governed by the thermal diffusion coefficient of intrinsic tissue and microvascular architecture rather than by bulk convective transport in vessels.

3.3 Characterization of external-surface convective heat transfer

The effective external convective heat-transfer coefficient, h, is a key design parameter that sets the boundary heat flux according to the equations (1) and (3), thereby governing the overall cooling efficiency, thermal non-uniformity, and the time required for the organ to reach near-equilibrium temperature^[48]. Figure 5 summarizes the sensitivity of predicted cooling to h for the liver, heart, and kidney. Figure 5a,c,e show that increasing h accelerates cooling most noticeably when moving from 1 (W/m²·K), while further increases to 100-500 (W/m²·K) produce comparatively smaller changes. This indicates diminishing returns once external convection is no longer the dominant resistance when h is over 100 (W/m²·K).

Display Full Size

Figure 5. Parametric analysis of external-surface convective heat-transfer coefficients for supercooled machine perfusion of the liver, heart, and kidney. (a, c, e) Simulated temporal evolutions of representative parenchymal temperature for the (a) liver, (c) heart, and (e) kidney, respectively, under different h; (b, d, f) Final time-point organ parenchymal temperature slices showing the spatial temperature distribution of the (b) liver, (d) heart and (f) kidney, respectively, under different h. Color legends indicate T (K).

The final-time temperature slices (Figure 5b,d,f) further illustrate that low h leads to higher overall tissue temperatures and stronger spatial non-uniformity near the external surface, consistent with limited heat removal at the boundary. As h increases, the organ exterior becomes more effectively coupled to the cold environment, reducing warm surface regions and producing a more uniformly cooled parenchyma at the final time point. The ΔT_max variation indicated that increasing h improved intra-organ temperature uniformity (Figure S5d,e,f). Collectively, these observations support that ensuring h is over 100 (W/m²·K) reduces temperature gradients of organ parenchyma and enables rapid cooling, thereby minimizing thermal stress damage to the greatest extent possible.

3.4 Parameterization of perfusate specific heat capacity

The thermophysical properties of the perfusate exert a direct and substantial influence on parenchymal temperature uniformity and cooling rate^[49]. Therefore, we systematically examine the impact of c_p on transient thermal response during supercooled MP. The perfusate c_p governs the thermal storage capacity of the circulating medium, thereby sustaining a stronger driving temperature gradient between the tissue and the perfusate throughout the vascular network.

Figure 6 summarizes the sensitivity of predicted cooling to c_p for the liver, heart, and kidney under identical inlet and environmental conditions. In the time histories shown in Figure 6a,c,e, increasing c_p consistently enhances cooling performance. Higher specific heat capacity reduces the temperature rise of the perfusate as it traverses the organ, thereby thereby maintaining a larger temperature difference available for heat extraction from the surrounding tissue. Consequently, organs perfused with higher c_p media achieve lower parenchymal temperatures more rapidly and approach the bath temperature sooner. Notably, the separation among the temperature curves diminishes at higher c_p values. This indicates that beyond a threshold (4,000 J/(kg·K)), further increases yield diminishing improvements in cooling rate because the process becomes increasingly constrained by internal thermal diffusion and architecture rather than by perfusate heat capacity alone.

Display Full Size

Figure 6. Parametric analysis of perfusate specific heat capacity for supercooled machine perfusion of the liver, heart, and kidney. (a, c, e) Simulated temporal evolutions of representative parenchymal temperature for the (a) liver, (c) heart, and (e) kidney, respectively, under different values of c_p; (b, d, f) Organ parenchymal temperature slices at the final time point showing the spatial temperature distributions of the (b) liver, (d) heart, and (f) kidney, respectively, under different c_p. Color legends indicate T (K).

Figure 6b,d,f show that the spatial temperature maps at the final time point under lower c_p conditions are associated with warmer residual regions and more pronounced internal gradients, consistent with faster perfusate warming rate and reduced local heat-removal capacity. Quantitative metric figure further show that increasing c_p reduced ΔT_max, enabling the perfusate to remove more heat from the organ and improving intra-organ temperature uniformity (Figure S5g,h,i). Collectively, these findings indicate that perfusate c_p is a critical perfusion parameter enhances cooling efficiency, improves thermal uniformity, and shortens the time to near-equilibrium conditions, while c_p over 4,000 J/(kg·K) produces limited gains due to intrinsic tissue transport limitations.

3.5 Effect of perfusate thermal conductivity

In addition to the perfusate heat capacity, its thermal conductivity k directly modulates the efficiency of heat exchange between the vascular lumen and surrounding parenchyma by altering the local conductive resistance within the flowing medium and at the tissue–perfusate interface^[50]. k becomes a key thermophysical property when convection is strong and the remaining limitation shifts toward conduction across small length scales and complex microvascular pathways^[51].

Figure 7 summarizes the sensitivity of predicted cooling to perfusate thermal conductivity for the liver, heart, and kidney under identical inlet and environmental conditions. The representative temperature time histories (Figure 7a,c,e) indicate that increasing k improves cooling performance, with the most pronounced benefit observed when k is increased from very low values (0.01-0.1 W/(m·K)) to more solution-like values (0.3-1 W/(m·K)), as shown in Figure 7b. Physically, a higher k enhances radial heat conduction within the perfusate and reduces the effective thermal resistance for transferring heat out of the perivascular tissue, enabling faster equilibration between the vascular domain and adjacent parenchyma. However, the separation among the curves progressively diminishes at higher k (≈1-2 W/(m·K)), indicating diminishing returns once interfacial/perivascular conductive resistance is no longer the dominant bottleneck, and the cooling dynamics become increasingly governed by intrinsic tissue thermal diffusion and microvascular architecture, an interpretation consistent with prior cryopreservation heat-transfer analyses emphasizing internal transport limits and thermo-mechanical consequences in large specimens.

Display Full Size

Figure 7. Parametric analysis of perfusate thermal conductivity for supercooled machine perfusion of the liver, heart, and kidney. (a, c, e) Simulated temporal evolutions of representative parenchymal temperature for the (a) liver, (c) heart, and (e) kidney, respectively, under different values of k; (b, d, f) Organ parenchymal temperature slices at the final time point showing the spatial temperature distributions of the (b) liver, (d) heart, and (f) kidney, respectively, under different k. Color legends indicate T (K).

The final-time temperature maps (Figure 7b,d,f) reinforce this mechanism: low k perfusates exhibit warmer residual regions and more pronounced intra-organ gradients, consistent with poorer conductive coupling between the cold perfusate and surrounding tissue. As k increases, these warm pockets contract, and the parenchymal field becomes more spatially uniform, reflecting more effective distribution of cooling capacity throughout the vascular network. This trend was further supported by the quantitative metric ΔT_max. Figure S5j,k,l shows that ΔT_max decreased with increasing k, indicating improved intra-organ temperature uniformity. Collectively, these results identify perfusate thermal conductivity as an actionable parameter for improving cooling efficiency and thermal uniformity, while also highlighting that pushing k beyond 1 W/(m·K) yields limited additional benefit because organ-scale cooling becomes constrained primarily by tissue-side transport rather than perfusate-side conduction.

3.6 Sensitivity analysis

We performed a one factor at a time sensitivity analysis of h and the perfusate thermophysical properties, including c_p and k, to identify the parameters that most strongly govern cooling during subcooled MP. The sensitivity analysis used the terminal parenchymal temperature and the cooling rate at the liver center as the primary response metrics. Each parameter was varied by ±20% about its baseline value, while all other model inputs were held constant.

Figure 8 shows that the influence of c_p on the cooling process in supercooled mechanical perfusion substantially exceeds that of h and k. Specifically, for the kidney, the changes in terminal temperature (> 0.01 K) and cooling rate (> 0.008 K min^-1) caused by variations in the heart’s c_p were substantially larger than the corresponding changes attributable to h and k (Figure 8c,d). Additionally, this conclusion is contingent on the baseline assumption h = 100 W/(m²·K). Therefore, adequate convection in the storage chamber remains essential for effective external heat removal.

Display Full Size

Figure 8. Sensitivity of organ cooling outcomes to convective heat transfer and perfusate thermophysical properties during supercooled machine perfusion. (a, b) Liver; (c, d) Heart; (e, f) Kidney. Panels (a, c, e) show changes in terminal parenchymal temperature, and panels (b, d, f) show changes in cooling rate; k, c_p and h parameters were individually perturbed by ±20% about their baseline values while all other inputs were held constant; (e, f) The change of k is too small to be clearly visible.

Collectively, when the external h at the organ surface exceeds 100 W/(m²·K), improving the thermophysical properties of the perfusion solution can provide substantial gains in cooling performance during supercooled MP. Accordingly, optimizing supercooled MP protocols in multiorgan preservation may enhance cooling performance by increasing the c_p of the perfusion solution.

4. Conclusions

This study presents a comprehensive heat-transfer analysis of supercooled perfusion for the liver, heart, and kidney using an anatomically informed porous-media modeling framework under experimentally relevant operating conditions. By integrating realistic perfusion boundary conditions with three-dimensional organ geometries derived from the BodyParts3D repository, the model successfully captured transient intra-organ temperature evolution during supercooled perfusion. Experimental validation using a self-developed machine perfusion platform demonstrated good agreement between simulated temperature trajectories and thermocouple measurements, with MAEs of 0.24 °C for the liver, 2.63 °C for the heart, and 0.4 °C for the kidney, supporting the reliability of the proposed modeling approach. The results revealed pronounced organ-dependent thermal responses to perfusate cooling, with larger organs, such as the liver, exhibiting slower temperature transients compared with the heart and kidney.

Effective cooling rates were found to decrease with increasing organ size and to increase with higher perfusion flow rates, highlighting the critical role of perfusion parameter selection in achieving efficient and uniform cooling. However, further increases in perfusion rate produced only minor changes in temperature uniformity. Accordingly, the suitable inlet conditions were set to a liver perfusion pressure of 30 mmHg, a heart inlet velocity of 1 mm/s, and a kidney perfusion pressure of 50 mmHg. These findings underscore the value of computational modeling as a quantitative tool for thermal design and analysis of supercooled perfusion protocols, particularly in multi-organ preservation scenarios where direct intra-organ temperature measurements are limited. Furthermore, sensitivity analysis further identified the c_p of perfusion solutions as the optimal parameter for supercooled MP protocols.

Despite its demonstrated predictive capability, the present model adopts a homogenized porous-media representation of tissue and simplified perfusion assumptions, which may not fully capture local heterogeneities in vascular architecture or perfusion distribution, especially in complex organs such as the heart. Additionally, unmodeled vessels and experimental uncertainties may contribute to minor deviations between simulated and measured temperatures. Moreover, because temperature measurements were obtained at only limited representative locations, the current validation primarily supports local temperature evolution and overall cooling trends rather than a fully resolved three-dimensional temperature field. Accordingly, the current framework is expected to provide more reliable temperature predictions for the kidney and liver under conditions broadly consistent with the present modeling assumptions and operating conditions. In contrast, its application to strongly anisotropic organs, such as the heart, where detailed vascular architecture may play a dominant thermal role, should be interpreted with caution.

In summary, this study establishes a validated, under the tested conditions computational framework for predicting intra-organ temperature evolution during supercooled MP and provides a theoretical foundation for thermal optimization of perfusion strategies aimed at improving preservation efficacy in multi-organ systems. The framework is reliable within the examined range of organ types, perfusion conditions, and detailed vascular architectures, while extrapolation beyond these conditions should be undertaken cautiously.

Future work should focus on incorporating more detailed vascular representations, organ-specific perfusion characteristics, and coupled thermo-fluid-biological effects to further enhance model fidelity. Systematic optimization of perfusion parameters across a wider operating space may also provide deeper insight into protocol design for large-scale organ preservation.

Supplementary materials

The supplementary material for this article is available at: Supplementary materials.

Authors contribution

Liu Z: Writing-original draft, investigation.

Wang X: Writing-original draft, data curation.

Zhang J: Investigation.

Yang S: Investigation, data curation.

Zhang H: Project administration.

Li Z, Wang Y, Yang F, Hou Y, Guo M: Data curation.

Rao W: Writing-review & editing, resources, methodology, funding acquisition.

Conflicts of interest

Wei Rao is an Editorial Board Member of Thermo-X. The other authors declare no conflicts of interest.

Ethical approval

All animal experiments and care were approved by the Animal Ethics Committee of the Medical Service Biotechnology Co., LTD (Approval No. MDSW-2025-031C).

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

Data supporting the findings of this study are available from the corresponding author upon reasonable request.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 52406117), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB1030000), Fund of State Key Laboratory of Cryogenic Science and Technology (Grant No. E3ASR201), Beijing United Fund (Grant No. L252063&L252064), and TIPC director fund (Grant No. E4AHR021).

Copyright

© The Author(s) 2026.