Intermolecular interactions play a role in the distribution and transport of charged contrast agents in a cartilage model

Summary

The transport and distribution of charged molecules in polyelectrolyte solutions are of both fundamental and practical importance. A practical example, which is the specific subject addressed in the present paper, is the transport and distribution of charged species into cartilage. The charged species could be a contrast agent or a drug molecule involved in diagnosis or treatment of the widespread degenerative disease osteoarthritis, which leads to degradation of articular cartilage. Associated scientific issues include the rate of transport and the equilibrium concentrations of the charged species in the cartilage and the synovial fluid. To address these questions, we present results from magnetic resonance micro-imaging experiments on a model system of articular cartilage. The experiments yield temporally and spatially resolved data on the transport of a negatively charged contrast agent (charge = -2), used in medical examinations of cartilage, into a polyelectrolyte solution, which is designed to capture the electrostatic interactions in cartilage. Also presented is a theoretical analysis of the transport where the relevant differential equations are solved using finite element techniques as well as treated with approximate analytical expressions. In the analysis, non-ideal effects are included in the treatment of the mobile species in the system. This is made possible by using results from previous Monte Carlo simulations. The results demonstrate the importance of taking non-idealities into account when data from measurements of transport of charged solutes in a system with fixed charges from biological polyelectrolytes are analyzed.

Keywords:

magnetic resonance imaging – cartilage – Mass diffusivity – Finite element analysis – Solute transport – Convection – Electrostatics

Introduction

The function of articular cartilage is to provide low friction and wear resistance, as well as to distribute load in load-bearing joints, for instance, in hips and knees [1]. The optimal functioning of the joint depends on the structure, composition as well as the integrity of the extracellular matrix in cartilage, which is mainly composed of water, collagen and proteoglycans [1, 2]. The load resistance property of cartilage is largely dependent on the proteoglycan content, originating from the glycosaminoglycan (GAG) side chains of the proteoglycan [3]. The GAG is highly negatively charged on account of its carboxyl and sulfate groups that are ionized under physiological conditions. These negative charges are fixed within the extracellular matrix and thus the expression fixed charge density (FCD) of the cartilage is used. Please note that the value of FCD is negative.

Osteoarthritis is a degenerative disease that affects the function of cartilage. The disease is a common cause of disability and it is therefore of importance to understand its progression and to develop efficient treatments. In early stages of the disease, a loss of GAG occurs [4] and one method used to monitor the resulting decrease in the FCD is delayed gadolinium-enhanced magnetic resonance imaging of cartilage (dGEMRIC) [5]. This method utilizes magnetic resonance imaging (MRI) together with an intravenous injection of a contrast agent. The contrast agent commonly used is Gd(DTPA)²⁻ which is negatively charged and contains the paramagnetic ion Gd³⁺, which has a concentration dependent influence on the spin-lattice relaxation time, T₁, of the ¹H of the water [6]. The dGEMRIC method is based on the fact that the Gd(DTPA)²⁻ complex will, due to its negative charge, distribute in an inverse relation to the GAG content of the cartilage. This means that where there is a high concentration of GAG, indicating a healthy cartilage, there will be a low concentration of contrast agent, and vice versa. Therefore, a measure of the concentration of the Gd(DTPA)²⁻, and consequently a proxy for the GAG content of cartilage, is obtained by measuring T₁ of the water. In passing, we note that the dGEMRIC method is used both in vitro and in vivo [5, 7–12]. The method has found most of its applications in in vitro studies and this will most likely continue to be the most important application.

dGEMRIC is dependent on the transport of Gd(DTPA)²⁻ from the site of the injection to the joint of interest, and, subsequently, on the diffusion controlled transport into cartilage, since cartilage lacks direct blood supply. It is therefore important to understand the physiochemical mechanisms that govern the transport and partitioning of Gd(DTPA)²⁻ in cartilage. The transport of the contrast agent will, for example, be affected by the chemical potential difference between the cartilage and the synovial fluid and, in this respect, the FCD of the cartilage plays an important role. Furthermore, the transport will be affected by the amount of collagen and glycosaminoglycan on account of steric effects [13–17]. This is an example of the general and important problem of diffusion of small molecules in a complex micro-heterogeneous environment, where crowding and excluded volume effects as well as electrostatic effects are important. The results presented here are therefore applicable to a wider class of related problems than just issues dealing with cartilage.

A previously developed magnetic resonance micro-imaging (μMRI) setup for the investigation of features related to the dGEMRIC method in a model system of cartilage was used to investigate the questions raised above under controlled conditions [18]. The μMRI setup gives spatially and temporally resolved information on the transport of Gd(DTPA)²⁻ in the model system, which is designed with the goal of investigating the effects of the electrostatic interactions in cartilage. To this end, the cartilage is represented by a polyelectrolyte solution and the synovial fluid is represented by a NaCl solution. A semipermeable membrane is used to separate the polyelectrolyte from the salt solution. The experimental results are compared with numerical solutions obtained using the finite element method (FEM) of the relevant equations describing the transport of ionic species in our system. Since previous studies have shown that ions in cartilage do not behave ideally [19–21], the influence of intermolecular interactions are included in the theoretical work presented in this paper. This is made possible by comparing the obtained data with previous Monte Carlo simulations on the model system [21]. Furthermore, the low concentration of Gd(DTPA)²⁻ makes it possible to describe the transport with an approximate analytic solution for a system without boundaries.

In summary, our research goal is to address the transport of charged molecules in complex micro-heterogeneous systems and to develop methods where intermolecular interactions are accounted for.

Model system

The model system is outlined in Fig 1A–1C. There is a polyelectrolyte solution at z < 0 and a salt solution at z > 0, separated by a semipermeable membrane at z = 0. The charged polyelectrolyte cannot pass through the membrane into the salt solution, which leads to a difference in electric potential on the two sides of the membrane, which in turn will be compensated for by a concentration difference of ions. The following scenario with respect to time, t, is investigated:

• At t < 0, there is equilibrium of Na⁺ and Cl⁻ ions between a salt solution and a solution containing a negatively charged polyelectrolyte.

• At t = 0, a small amount of Gd(DTPA)²⁻ (see Fig 1D for the chemical structure of the complex) is added to the salt solution. The concentration of Gd(DTPA)²⁻ is much lower than the concentration of the other charged species in the system.

• At t > 0, the Gd(DTPA)²⁻ will diffuse into the polyelectrolyte solution.

The concentration difference of Gd(DTPA)²⁻ between the solutions at t > 0 will lead to transport of the Gd(DTPA)²⁻ into the polyelectrolyte solution. Initially, the difference is large and leads to fast transport by mutual diffusion from the salt solution into the polyelectrolyte solution (a flow of Na⁺ and Cl⁻ ions will also occur to maintain the electro-neutrality in the system, which is minor compared to the total concentration of Na⁺ and Cl⁻ ions). As equilibrium is approached, the rate of transport decreases. Since the transport is only driven by diffusion, it will take a significant length of time to reach equilibrium. Recall that the root mean square displacement is proportional to the square root of time for diffusive motion.

Theoretical background

The basic equations used in the interpretation of the data are presented in this section. The full derivation of the relations used in the Results and Discussion section is presented in the accompanying S1 and S2 Appendices.

The diffusive flux of an ionic species i, J_i, occuring during step 3 above, can be described by the generalised Fick’s first law (in one dimension):

The electric potential, ϕ, is related to the charge distribution through Poisson’s equation:

To be able to solve Eqs 3 and 4, μ i corr also needs to be known. It can be estimated by using results from previously made Monte Carlo simulations performed on a system developed to capture the electrostatic interactions in the cartilage/synovial fluid system at equilibrium [21]. In the simulations, the following holds for the two separate solutions, if the electric potential in the salt solution is set as the reference:

Material and methods

Materials

Sodium carboxymethylcellulose (CMC) (average molecular weight = 90 kDa, 99.5% purity, Sigma Aldrich; molecular structure is given in Fig 1E) was used without further purification. The degree of substitution of the used batch of the polyelectrolyte has previously been determined to 0.9 [18]. Magnevist^® (0.5 mmol gadolinium/mL, Bayer Healthcare) was used as the Gd(DTPA)²⁻ source. The molecular structure of Gd(DTPA)²⁻ is given in Fig 1D. NaCl (99.5% purity, Merck) and NaOH (analytical grade) were used as received. All solutions were prepared with MilliQ purified H₂O. Regenerated cellulose acetate membranes with a cutoff of 5 kDa were obtained from Millipore. This cutoff was previously shown not to affect the transfer rate of Gd(DTPA)²⁻ significantly and at the same time only allow for a negligible amount of CMC to pass through the membrane [18]. For dialysis of the polyelectrolyte solutions, Spectra/Por^® Biotech Cellulose Ester (CE) Dialysis membrane (dialysis cutoff of 100-500 Da) was used.

Sample preparation

The pH of the salt solution (150 mM NaCl) was adjusted to between 8 and 9 by means of drop-wise addition of 1 M NaOH. This pH value is slightly higher than the generally quoted value of synovial fluid of around 7.5. In this pH range, effectively all of the carboxylic groups on CMC are charged.

Magnevist^® was diluted by the salt solution to a Gd(DTPA)²⁻ concentration of 0.36 mM. Before use, the pH was adjusted to the same value as in the salt solution.

When preparing the polyelectrolyte solutions, the CMC powder was dissolved in the salt solution (150 mM NaCl). Before use, the polyelectrolyte solutions (7 mL) were dialyzed for 24 h against 1.5 L 150 mM NaCl solution. A control experiment showed that this time was enough to reach ionic equilibrium between the polyelectrolyte solution and the salt solution (for more information, see S3 Appendix). The salt solution (the dialyzate) was gently stirred during the dialysis and the pH of the dialyzate was kept between 8 and 9 by periodically making drop-wise addition of 1 M NaOH. The carbon content in the polyelectrolyte solutions was determined by total organic carbon (TOC) analysis after the dialysis. The FCD was then calculated using as input the degree of substitution of the polyelectrolyte. To obtain further information with regard to the ionic content, the sodium content was determined in the solutions by inductively coupled plasma atomic emission spectroscopy (ICP-AES) analysis. Subsequently, the Cl⁻ concentration could be determined by applying electroneutrality conditions in the solutions. The densities of the solutions were determined by a precision density meter (DMA5000, Anton-Paar). The polyelectrolyte concentration of the polyelectrolyte solutions was estimated to 1.7, 2.2 and 2.5 wt%, respectively.

A small piece of the regenerated cellulose acetate membrane was punched out to fit into the sample holder (see Fig 1B and 1C). Before use, the membrane was, according to instructions from the manufacturer, put in 150 mM NaCl solution (with the glossy side up) to wash out glycerol and to wet the pores. The solution was changed three times during 1 h and thereafter, the membrane was left in the solution for at least 24 h before use.

Experiments

The bottom section of the sample holder (constructed from polyether ether ketone (PEEK) plastic and a 10 mm NMR tube; see Fig 1B and 1C) was filled with ∼160 μL CMC solution (added by weight). Thereafter, the membrane was carefully placed on top of the polyelectrolyte solution with the glossy side down. Subsequently, the top part of the sample holder (including the NMR tube) was assembled with the bottom part and 205 μL of the 150 mM NaCl solution was added on top of the membrane.

2.1 mL of the 0.36 mM Gd(DTPA)² solution was added to the top compartment about 30 minutes after the addition of the salt solution (see above). The first T₁ measurement (for description of T₁ measurements, see below) was started within a few minutes after the addition of Gd(DTPA)²⁻. Due to rapid initial changes of the Gd(DTPA)²⁻ concentration, the initial T₁ measurement (for each sample) gave unphysical concentration profiles in the salt solution and these data points were not used in the subsequent analysis. The sample holders were kept in the spectrometer for the initial 22 h and T₁ measurements were performed continuously. Subsequently, the sample holders were removed from the spectrometer. For acquisition of subsequent data, the sample holder was re-inserted into the spectrometer for each measurement. Data at 3 different positions were collected in each measurement (one set for the whole polyelectrolyte solution, one set simultaneously showing the polyelectrolyte and salt solutions, and one set for the salt solution close to the membrane). This measurement procedure was carried out once a day. Between the measurements, the sample holders were stored at 298 K, protected from light. The samples were not stirred during the experiments and, in order to minimize agitation, the sample holders were carefully handled between measurements. Three different concentrations of CMC were studied and for each concentration, two independent experiments were performed.

μMRI measurements

Spatially resolved T₁ measurements were recorded at 298 K on an 11.7 T Bruker Avance II 500 spectrometer equipped with a Bruker MIC-5 microimaging probe having a maximum gradient strength of 3 T m^-1 and a 10 mm saddle-coil radio-frequency insert. The measurements were performed using a pulse sequence containing an inversion recovery block followed by a rapid acquisition with relaxation enhancement (RARE) block [22] for the read out of the images. The recovery time delay, τ, was incremented logarithmically in 32 steps from 0.0017 s to 15 s, while the images were taken up with an 11.1 × 4.8 mm field-of-view (z × y), 128 × 48 acquisition matrix size and a 0.2 mm slice thickness (x). This resulted in a measurement time of about 10 minutes per image.

Data processing

The images were reconstructed to a 128 × 48 matrix size, which gives a spatial resolution of 87 × 100 μm, and subjected to 0.2 mm Gaussian smoothing. For each voxel, the T₁ values were obtained by fitting the image intensity, I, to following equation:

The obtained T₁ values were then converted to the concentration of Gd(DTPA)²⁻, c G . 2 −, using the expression:

The three different images that were taken at the same time point (see above) were visually merged to obtain a larger field-of-view of the system.

Assuming that there is no convection in the polyelectrolyte solution, the diffusion coefficients of Gd(DTPA)²⁻ in the polyelectrolyte solutions were obtained from fitting the relevant equations to the concentration profiles in the polyelectrolyte solution (see S1 Appendix for explanation of the procedure). The data points used in the fit were obtained in a time interval where the influence of the membrane is low and the Gd(DTPA)²⁻ has not yet reached the bottom of the polyelectrolyte compartment, i.e. 2.5-6 h.

All data processing was performed with in-house written Matlab^® code.

Numerical methods

Numerical calculations were carried out using COMSOL Multiphysics^® 4.4, a program that solves partial differential equations using the finite element method (FEM). The FEM simulations were carried out for a one-dimensional geometry with the same dimensions as in Fig 1C, except that the salt solution was truncated at z = d_salt to account for thermally driven convection in the salt solution. At this distance, the concentration in the salt solution in the FEM simulations was set equal to the experimental value at the same position. The distance d_salt was varied between 1 to 5 mm for the different experiments, depending on whether the contribution from convection was high or low, respectively. This was determined by comparing the experimental concentration profiles in the salt solution with those obtained from the FEM simulations without convection and choosing d_salt at a value where the experimental and the calculated curves in the salt solution are equal. Moreover, we assume homogenous properties in each solution, such that values of the diffusion coefficients are constant in each solution and the chemical potential correction factors are, except close to the interface, also constant in each solution.

Data from all experiments (two sets of data for each polyelectrolyte concentration) were used in the simulations. First, the time-independent Eq 4 and Poisson’s equation (Eq 6) were solved for the Na⁺ and Cl⁻ ions as well as for the electric potential. In order to avoid a singularity in the change in chemical potential at z = 0, we introduced a smooth change in μ i corr around z = 0 according to:

The boundary conditions used are summarized in Table 2. The initial concentrations were c Na + , s = 150mM, c Cl − , s = 150mM, c G . 2 − , s = 0.33mM and c G . 2 − , p = 0mM. The temperature was 298.15 K, ϵ_r = 78, and the diffusion coefficients for Na⁺ and Cl⁻ were assigned values of 1.33⋅10^-9 m²/s and 2.03⋅10^-9 m²/s, respectively [24]. It should be noted that the obtained results are independent of the chosen diffusion coefficients, since we are assuming a situation of steady state for the Na⁺ and Cl⁻ ions. The value for the diffusion coefficient for Gd(DTPA)²⁻ used in the salt solution was obtained from a study of molecules similar to Gd(DTPA)²⁻ and was taken to be 4.5⋅10^-10 m²/s [25]. The diffusion coefficient for Gd(DTPA)²⁻ in the membrane was set to 1.5⋅10^-10 m²/s in all FEM simulations, chosen to fit with the experimental data at early measurement times. As noted above, the diffusion coefficients in the polyelectrolyte solution were obtained from analysis of the experimental data (see also below and in the S1 Appendix). Values for Δ μ i corr were determined from the equilibrium concentrations together with equations describing equilibrium concentrations in the framework of a Donnan equilibrium. However, values for Δ μ i corr cannot be determined independently for Na⁺ and Cl⁻ as there are two independent relations and three unknowns (Δ μ Na + corr, Δ μ Cl − corr and ΔμG.2−corr, see Equations D-F in S2 Appendix). To arrive at individual values of Δ μ i corr needed in the analysis of Gd(DTPA)²⁻, we start by setting the Δ μ Na + corr to values of Δ μ Na + ex obtained from Monte Carlo simulations [21]. The Δ μ Cl − corr values are subsequently obtained from the equilibrium concentrations of Na⁺ and Cl⁻ (see Table 1) and Equations D and E in S2 Appendix. In principle, values for Δ μ G . 2 − corr can then be determined from Equation F in S2 Appendix. However, the equilibrium concentrations of Gd(DTPA)²⁻ are less accurately known, on account of the long time required to reach equilibrium. To improve the accuracy in the obtained value, we have made an initial estimate from Equation F in S2 Appendix and then allowed the equilibrium value of Gd(DTPA)²⁻ to vary in the FEM simulations so that the simulated concentration profiles best fit the experimentally determined ones. The used values are collected in Table 1.

To investigate the validity of the approximate equations that are derived in the Results and Discussion section, FEM simulations were also carried out with same dimensions as in Fig 1C (not truncated at a distance d_salt as described above) and values taken from Monte Carlo simulations (Ref. [21] and unpublished results, see S3 Fig).

Results and discussion

The outline of this section is as follows. We start by presenting and discussing the experimental data obtained, and then proceed to numerically solve the relevant equations describing the data using FEM. Subsequently, we develop approximate solutions that can be used to illustrate how different parameters affect the distribution of ions in the system. Using examples, we demonstrate under which conditions the approximate relations are valid. Finally, we use the approximate relations to determine the diffusion coefficients of the contrast agent in the polyelectrolyte solution. In the interest of keeping the paper short, we delegate derivations to the S1 and S2 Appendices, to which we refer when needed.

Experimental data

The main experimental results of this study are concentration profiles for Gd(DTPA)²⁻ as a function of time for three different values of FCD. An example is given in Fig 2, which shows the experimentally determined concentrations of Gd(DTPA)²⁻ in a system with the polyelectrolyte solution having a FCD of -108 mM, at different times after addition of Gd(DTPA)²⁻ to the salt solution (results for the duplicate experiment with FCD = -108 mM and the other values of FCD are shown in the S1 Fig). After the addition, Gd(DTPA)²⁻ diffuses from the salt solution to the polyelectrolyte solution. This is a slow process as can be observed in Fig 2. The times needed for the average concentration in the polyelectrolyte solution to reach 50% and 90% of the equilibrium concentration in the 10 mm long polyelectrolyte solution are 40 h and 230 h, respectively. However, there is, as expected, a rapid development of a discontinuity in the concentration at z = 0, although it appears to be gradual in Fig 2 due to the finite resolution of the μMRI method and artifacts due to the presence of the membrane (e.g. uncertain values of relaxivities and T_1pre). The FEM simulations make it possible to study the discontinuity closer (see below).

In Fig 3A, we show the concentration of Gd(DTPA)²⁻ in the polyelectrolyte solution at a position 4 mm from the membrane as a function of time for a FCD of -108 mM. The concentration of Gd(DTPA)²⁻ will initially increase rapidly in the polyelectrolyte solution, but the magnitude of the increase will gradually decrease with time. The symbols present two sets of data recorded in repeated experiments with the same FCD for the polyelectrolyte solution. The difference in the curves is attributed to a different degree of temperature gradient driven convection in the salt solution at the beginning of the experiments. The contribution from convection can be estimated by analyzing the concentration profiles in the salt solution (see Fig 2A and S1 Fig). Furthermore, convection evens out concentration gradients in the salt solution and therefore, the time to reach equilibrium becomes shorter. On account of the geometry (small radius and height) and the position of the sample holder in the spectrometer (where temperature gradients are small), and the fact that the polyelectrolyte is rather viscous, the transport in the polyelectrolyte solution is well described by pure diffusion.

In the following sections, we will use relevant theoretical expressions in combination with FEM simulations to investigate and rationalize the behavior of the Gd(DTPA)²⁻ with regard to its distribution and transport. The analysis yields physical properties relevant for the system. We note that the approach can be used to predict and understand the transport of ions in general in related systems.

Finite element method simulations

By solving Eqs 3, 4 and 6 numerically, using FEM techniques (as described above in Material and methods), we are able to predict the experimental results in Fig 2. Using as input the data in Table 1 (and values for the diffusion coefficients for Na⁺ and Cl⁻ given above) and the boundary conditions in Table 2, the results are presented as solid lines in Fig 3A and 3B (for additional comparison between experiments and FEM simulations, see S2 Fig). As can be seen, there is a good agreement with the experimental results at all times, and the effects due to convection in the salt solutions are well accounted for. We stress that all relevant parameters have been determined in independent experiments, except for the values of the electric potential and the chemical potential correction factor for the mobile ions, which are in good agreement with the results from independent Monte Carlo simulations (Ref. [21] and unpublished results, see S3 Fig).

One of the key motivations for carrying out this study was to investigate to what extent effects due to non-ideal conditions are important in the transport and equilibrium concentrations of Gd(DTPA)²⁻ in cartilage and models of cartilage. This is one important issue in the dGEMRIC method. To this end, we present in Fig 3A and 3C results of FEM simulations where the difference in the chemical potential correction factors between the two solutions of all species are set to zero. It is clear that the predicted equilibrium concentrations deviate considerably, in agreement with earlier work [5, 18, 21].

The FEM simulations enable us to study the discontinuity in the concentration at z = 0. In Fig 3D, the step in the Gd(DTPA)²⁻ concentration, Δ c Gd ( DTPA ) 2 − (see also Fig 4), is plotted against time. There is a rapid development of a discontinuity in the concentration at z = 0 and the magnitude of the step increases with time. The double experiments for each value of FCD, show that the magnitude of the step size is reproducible. The deviations are due to different degrees of convection in the salt solution. While the steps in concentration vary with time, the ratios between the concentrations of Gd(DTPA)²⁻ in the polyelectrolyte solution and salt solution at z = 0 are independent of time, with values of 0.82, 0.78 and 0.74 at FCD values of -73 mM, -92 mM and -108 mM, respectively.

Approximate relations for concentration profiles

Motivated by the need to illustrate the influence of the relevant parameters and to carry out rapid and straightforward calculations of the diffusion and concentration profiles, we here develop approximate expressions describing the transport and concentration profiles of Gd(DTPA)²⁻ in our model system. Our point of departure is to consider an interval 0⁻ ≤ z ≤ 0⁺ (see Fig 4), the size of which is such that we can assume that quasi steady-state conditions apply. This is justified by the fact that the gradient in the concentration at z = 0 adjusts much more rapidly than the variation in concentration over macroscopic length scales outside 0⁻ ≤ z ≤ 0⁺. In the current experiments, this interval is around 10 nm which is orders of magnitude smaller than the dimensions of the macroscopic container (see Fig 4). Solving Eq 3 with J_i = 0 gives (the reader is referred to S1 Appendix for details):

We next carry out an approximate analysis of the transport of Gd(DTPA)²⁻. To this end, we solve Eq 17 on both sides of the salt/polyelectrolyte solution interface at z = 0

Approximate analytic expressions vs finite element simulations

To illustrate under which conditions the approximate relations above are valid, we compare them with FEM solutions for the exact (within the model) equations. Fig 5 shows the concentration profiles using the analytical expressions in Eqs 18 and 19 and Equation W in S1 Appendix and how it compares to the simulated values at different times for a low FCD = -50 mM (Fig 5A), for the highest FCD used in this study, i.e. -108 mM (see Fig 5B), and for a higher FCD = -150 mM (see Fig 5C). It can be seen that the approximate analytic expressions reproduce the concentration profiles for values of FCD of -50 mM and -108 mM, respectively, but less so for FCD = -150 mM on account of the linearization of Equation B in S1 Appendix. We note that for short times (less than 2 h), it is non-trivial to include the effect of the membrane on the distribution of ions in the different compartments (see S1 Appendix for further information). The concentration profile in the polyelectrolyte solution starts to deviate considerably from the analytical expressions when the Gd(DTPA)²⁻ ions have reached the wall of the container containing the polyelectrolyte solution (not shown).

Diffusion coefficient of contrast agent in the polyelectrolyte solution

From Eq 18, it is possible to obtain values of the diffusion coefficient of Gd(DTPA)²⁻ in the polyelectrolyte solution from the experimentally determined concentration profiles. As outlined in S1 Appendix, one has to take into consideration the fact that Eq 18 applies to a system with infinite boundaries and that there is convection in the salt solution. When this is accounted for, Eq 18 is modified to Equation Z in S1 Appendix. A fit of Equation Z in S1 Appendix to one dataset is given in S1 Appendix, and the obtained diffusion coefficients are presented in Table 1. We note that the values follow the expected dependence on the magnitude of FCD. Thus, the values are 91, 87 and 82% of the value in the salt solution for FCD values of -73, -92 and -108 mM, respectively. As noted in the Introduction, steric effects may affect the diffusion coefficient of the contrast agent, which raises the question whether such specific effects differ between our model system and real cartilage. It is a general observation that electrostatic effects dominate diffusion of charged species in charged colloidal systems of different kinds. Moreover, results in Refs. [13–17] indicate that steric effects are relatively minor in real cartilage. Taken together these observations does not prove, but support the notion that the diffusion data for the contrast agents in our model system has bearing to the situation in real cartilage. Finally, we note that the problem of diffusion of charged species through charged colloidal systems is a general problem in a multitude of applications, and we feel that our suggested approach of obtaining such data is one of the main results of this work.

Conclusions

We have combined magnetic resonance imaging results with simulations on a previously developed model system with the aim of studying issues related to early diagnosis and in vitro and in vivo clinical studies of osteoarthritis. The flow of a charged contrast agent from a salt solution into a polyelectrolyte solution has experimentally been determined with spatial and temporal resolutions of 90 μm and 600 seconds, respectively. The relevant equations describing the flow have been solved using Finite Element Methods. The results clearly show that effects due to non-ideality must be taken into account when analyzing flow and equilibrium distributions of divalent species, such as contrast agents, in complex polyelectrolyte systems, such as cartilage. Approximate expressions for the flow have been derived and compared with numerical solutions of full equations in order to assess when the approximations are reasonable. Finally, the derived approximate expressions make it possible to determine the diffusion coefficient for the contrast agent in a polyelectrolyte solution from the experimental data.

Supporting information

S1 Appendix [pdf]
Diffusion-driven transport of ions in an infinite, one-dimensional Donnan system.

S2 Appendix [pdf]
Donnan equilibrium.

S3 Appendix [pdf]
Determination of dialysis time.

S1 Fig [pdf]
Additional experimental concentration profiles.

S2 Fig [pdf]
Additional concentration profiles from experiments and FEM simulations.

S3 Fig [pdf]
Obtained values of from Monte Carlo simulations.