Computational modeling of diffusive transport

Lately I’ve been working on building computational models which can be used to make predictions about the behavior of a diffusion experiment, to allow for the collection of large pools of data without the difficulty and time consumption associated with actually performing an experiment. This could prove to be a valuable tool for the continued development of our hemodialysis device, which has demonstrated excellent urea clearance in animal models but which may require some tuning to achieve the sharpest possible size cutoff, retaining BSA while clearing beta-2-microglobulin (or cytochrome C, as the case may be.)

I’ve implemented two different models which approximate the experimental data I’ve collected over the previous several months. The first is a finite element model in COMSOL Multiphysics, which includes a three-dimensional CAD model of the single-well CytoVu device and simulates the diffusion experiment given the diffusion coefficient of the relevant molecule(s) both in the bulk and in the membrane. This model is the more accurate of the two, but requires significant computation time (a single simulation takes on the order of fifteen minutes on my MacBook to reach completion.) The second model is purely analytical, and is lifted directly from Jess Snyder’s Ph.D. thesis. This model is one-dimensional and assumes equal well height for the filtrate and retentate, and as such its relatively low accuracy compared to the COMSOL model is unsurprising; however, this model is suitable for quickly acquiring large datasets and is highly predictive on small time scales.

The CAD model used in COMSOL is depicted below, in the form of an example graphic result from a simulation based on the diffusion coefficient of cytochrome C. The diffusion coefficient values used in this model can be acquired from another of Jess’s works, her “Sieving Programs” MATLAB package (specifically, finding_dmem.m,) which uses the Stokes-Einstein relationship to estimate the diffusion coefficient in water given the hydrodynamic radius of the molecule. It also derives a diffusion coefficient value for the molecule within the membrane from a modified form of the Renkin equation; however, as it turns out, in most cases the diffusion coefficient value inputted for the membrane has almost no effect whatsoever on the overall mass transport due to the relatively tiny thickness of the membrane compared to the total thickness of the system. Of course, experimental or literature values can also be used.

The results of the computational models for cytochrome C are shown on the graph below, along with experimental data for comparison. The takeaway from this is that the models are excellent for shorter time scales, but should be used with more caution as the system moves further along towards equilibrium. In this plot, “analytical” refers to Jess’ model, while “simulation” refers to COMSOL. The models appear to be similarly predictive for other molecules (urea and BSA; I have less experimental data with which to make comparisons for these molecules, and will likely collect more eventually.)

I’m toying with the idea of applying an empirical correction to one or both of these models — it may be as simple as subtracting a low-magnitude exponential or error function curve, although more likely there’s some nuance to consider. If I want to go down that road, I’d definitely want to collect more experimental data in order to measure the difference between the models’ predictions and reality for much longer time scales, so it may be an unnecessary time sink.

For now, these two models together provide a toolbox for making reasonable predictions of the behavior of diffusion. The obvious next step for application to hemodialysis is the addition of convective elements to the models. My current work concerns making the necessary modifications to these models to approximate a flow-counter flow dialysis device, as well as possibly non-dimensionalizing them in order to build a comprehensive model of the dialysis process. This could then be optimized along various parameters to determine the ideal characteristics of a real device, as well as used to preclude future bench-top and animal experiments which can be shown to be likely to fail using modeling.