Extracellular Space

Upon the identification of all no-production metabolites in the model the next step involves filling these gaps using minimally the three mechanisms described earlier. We first explore whether reaction directionality reversal and/or addition of reactions from Metacyc [18 (link)] absent from the original model links the problem metabolite with the present substrates. This is accomplished by using a database of candidate reactions consisting of (i) all reactions in the original model with their directionalities reversed and (ii) reactions from a curated version of the MetaCyc database including allowable transport mechanism entries between compartments (in the case of multi-compartment models). It should be noted here that all the reactions in the MetaCyc database are treated as reversible in the model. These reactions define set Database comprised of candidate reactions while set Model is composed of the original genome-scale model reactions. It should be noted that if none of the above two/three mechanisms is capable of connecting the cytosolic no-production metabolite in single/multi-compartment models then an uptake reaction is arbitrarily added to the model to restore connectivity. However, if a non-cytosolic metabolite (in the case of multi-compartment models) present in an inner compartment cannot be fixed by any of the above mechanisms it is flagged as unfixable given the employed mechanisms.
In addition to the binary variable w_ij defined previously, the proposed (GapFill) formulation relies on the binary variables y_j defined as follows:
yj={1if reaction j from the external database is added to the parent network0otherwise MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWG5bqEdaWgaaWcbaGaemOAaOgabeaakiabg2da9maaceqabaqbaeaabiGaaaqaaiabigdaXaqaaiabbMgaPjabbAgaMjabbccaGiabbkhaYjabbwgaLjabbggaHjabbogaJjabbsha0jabbMgaPjabb+gaVjabb6gaUjabbccaGiabdQgaQjabbccaGiabbAgaMjabbkhaYjabb+gaVjabb2gaTjabbccaGiabbsha0jabbIgaOjabbwgaLjabbccaGiabbwgaLjabbIha4jabbsha0jabbwgaLjabbkhaYjabb6gaUjabbggaHjabbYgaSjabbccaGiabbsgaKjabbggaHjabbsha0jabbggaHjabbkgaIjabbggaHjabbohaZjabbwgaLjabbccaGiabbMgaPjabbohaZjabbccaGiabbggaHjabbsgaKjabbsgaKjabbwgaLjabbsgaKjabbccaGiabbsha0jabb+gaVjabbccaGiabbsha0jabbIgaOjabbwgaLjabbccaGiabbchaWjabbggaHjabbkhaYjabbwgaLjabb6gaUjabbsha0jabbccaGiabb6gaUjabbwgaLjabbsha0jabbEha3jabb+gaVjabbkhaYjabbUgaRbqaaiabicdaWaqaaiabb+gaVjabbsha0jabbIgaOjabbwgaLjabbkhaYjabbEha3jabbMgaPjabbohaZjabbwgaLbaaaiaawUhaaaaa@99F0@
For the case of single compartment models, the task of identifying the minimal set of additional reactions that enable the production of a no-production metabolite i* is posed as the following mixed integer linear programming problem (GapFill).
Minimize∑j∈Databaseyj(GapFill) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaafaqabeqacaaabaGaemyta0KaemyAaKMaemOBa4MaemyAaKMaemyBa0MaemyAaKMaemOEaONaemyzau2aaabuaeaacqWG5bqEdaWgaaWcbaGaemOAaOgabeaaaeaacqWGQbGAcqGHiiIZcqWGebarcqWGHbqycqWG0baDcqWGHbqycqWGIbGycqWGHbqycqWGZbWCcqWGLbqzaeqaniabggHiLdaakeaacqGGOaakcqqGhbWrcqqGHbqycqqGWbaCcqqGgbGrcqqGPbqAcqqGSbaBcqqGSbaBcqGGPaqkaaaaaa@549F@
s.t


∑jwi∗j≥1 MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaadaaeqbqaaiabdEha3naaBaaaleaacqWGPbqAcqGHxiIkcqWGQbGAaeqaaOGaeyyzImRaeGymaedaleaacqWGQbGAaeqaniabggHiLdaaaa@3836@
∑jSijvj≥0∀i∈N MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaafaqabeqacaaabaWaaabuaeaacqWGtbWudaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabdAha2naaBaaaleaacqWGQbGAaeqaaOGaeyyzImRaeGimaadaleaacqWGQbGAaeqaniabggHiLdaakeaacqGHaiIicqWGPbqAcqGHiiIZcqWGobGtaaaaaa@3EF0@





In (GapFill), the objective function (13) minimizes the number of added reactions from the Database so as to restore flow through metabolite i*. Constraints (14) and (15) are identical to (6) and (7). Constraint (16) ensures that these additions are subject to a minimum of δ units for the no-production metabolite i* being produced. Constraint set (17), as in (GapFind), allows for the free drain of all cytosolic metabolites while bounds on reactions present in the Model are imposed by constraint set (18). Constraint set (19) ensures that only those reactions from the Database that have non zero flow are added to the model. This formulation restores flow through no-production metabolites in single compartment models. For multi-compartment models, the (GapFill) formulation is modified. First gaps in the cytosol are filled using the mechanisms described earlier for single compartment models. Specifically, the (GapFill) formulation is modified by replacing constraint (17) with constraints (11) and (12) reflecting the fact that no net production term can be imposed for metabolites present with compartments incapable of communicating directly with the extracellular space. The solution of formulation (GapFill) once for each no-production metabolite i* identifies one mechanism at a time for resolving connectivity problems in the model. It should be noted that through the use of integer cuts [29 ] multiple hypotheses can be generated to resolve these connectivity problems. In this study, we evaluate the merit of generated hypotheses and subsequently choose the most probable one using the following three criteria sequentially a) The added hypotheses should not have cycles: since the MetaCyc database consists of multiple copies of the same reaction (which are present in different organisms), there is a proclivity to fix metabolites by adding two copies of the same reaction in opposite directions (since all reactions in the MetaCyc database are considered reversible) thereby forming a cycle, b) We choose the hypotheses which enables production of the problem metabolite with the least number of modifications and c) We choose a hypotheses that has higher probability of being accurate based on our validation metrics (e.g., if two reactions are added, we choose the one with the better blast score). Note that a GAMS implementation of (GapFill) is available as an additional file [see Additional file 4].

Conductivity tensor images were reconstructed using an MRCI toolbox which is available at http://iirc.khu.ac.kr/toolbox.html (Sajib et al., 2017 (link)). The raw data was extracted from the k-space of the MR spectrometer. To minimize geometrical mismatches, the B1 phase maps and DWI were registered with the anatomical T2-weighted images after denoising and bias corrections. The B1 phase map, which is the spatial sensitivity distribution of the applied RF coil measured via MRI, is used to obtain σ_H (Katscher et al., 2009 (link)). From the MEMS data obtained after imaging experiment, the multiple echoes were combined to achieve a higher signal-to-noise ratio (SNR) using a weighting factor. The optimized phase maps were used to reconstruct σ_H (Gurler and Ider, 2017 (link)). The multi-b-value DWI data were corrected for eddy-current effects and geometrical distortions. The averaged images at b = 0 were linearly co-registered to the magnitude images of MEMS data, and the affine transformation matrix was used to non-linearly co-register the multi-b-value DWI (Smith, 2002 (link)). The conductivity of extracellular space can be defined as the product of ion concentration and mobility of charged particles. The following CTI formula was used for all conductivity tensor image reconstructions (Sajib et al., 2018 (link)): ${\sigma }_L=\alpha {\sigma }_e=\alpha \overline{c_e}{\mu }_e$ where σ_L is the low-frequency conductivity; α is the extracellular volume fraction; σ_e is the conductivity of extracellular space; $\overline{c_e}$ is the ion concentration; ${\mu }_e$ is the ion mobility. The apparent extracellular ion concentration $\overline{c_e}$ can be estimated as suggested by Sajib et al, (2018) (link). $\overline{c_e}=\frac{{\sigma }_H}{\alpha d_e^w+\left(1-\alpha \right)d_i^w\beta }$ where β is the ion concentration ratio of the intracellular and extracellular spaces; d_e^w and d_i^w are the extracellular and intracellular water diffusion coefficients, respectively. Since the ion concentration inside and outside of the giant vesicles are almost similar (β = 1), the low-frequency isotropic conductivity σ_L can be expressed as follows: ${\sigma }_L=\frac{\alpha {\sigma }_H}{\alpha d_e^w+\left(1-\alpha \right)d_i^w\beta }{d}_{\mathrm{e}}^w$
The details of the conductivity tensor reconstruction procedures followed those outlined in the works of Katoch et al, (2019) (link).

The initial conditions for the ion concentrations in the intracellular and extracellular domains are specified in Table 2. In the membrane (outside ion channels), all ion concentrations are set to zero. Furthermore, in the ion channels, the concentration of the ion species that is able to move through the channel is initially set up to vary linearly from the intracellular to the extracellular part of the channel, and the remaining ion concentrations are set to zero.
The background charge density, ρ₀, is set up such that the entire domain is electroneutral at t = 0. This means that ρ₀ = 0 everywhere except in the ion channels, where it is determined by the initial conditions for the ion species able to move through the channel.
In S1 Appendix, we show the simulation results for a different choice of ρ₀ than the default linear profile. More specifically, we consider initial conditions in the channels specified by a jump from the intracellular to the extracellular concentration in the center of the channel. The main results concerning the magnitude of the changes of the extracellular potential between the cells seem to be quite similar for this other choice of ρ₀ in the channels. However, some differences between the solutions are observed (see S1 Appendix).
In S1 Appendix, we also show the results of simulations with a different choice of initial conditions for the intracellular concentrations. More specifically, the intracellular Na⁺ concentration in cardiomyocytes is known to vary between species, with a concentration of about 4–8 mM for most mammalian species, but 10–15 mM in rat and mouse [31 (link)–34 (link)]. Thus, the default value of the intracellular Na⁺ concentration reported in Table 2 (12 mM) is mostly representative of rat or mouse cardiomyocytes, and a value of about 8 mM would be more realistic for healthy human cardiomyocytes [32 (link)]. In S1 Appendix, we compare the results of simulations using either 12 mM or 8 mM as initial conditions for the intracellular Na⁺ concentration. In addition, the intracellular Cl⁻ concentration is set to either 137.0002 mM or 133.0002 mM, respectively, in order to maintain electroneutrality at t = 0. In S1 Appendix, we observe that the two choices of initial conditions give very similar results.
The default geometry used in the simulations is specified in Fig 1 and Table 1. In the purely intracellular and extracellular spaces, the relative permittivity, ε_r, is set to ε₁ = 80. In the membrane and the ion channels, ε_r is set to ε_m = 2. Furthermore, in the purely intracellular and extracellular domains, Ω_i and Ω_e, the diffusion coefficients for the ions are as specified in Table 3. In the membrane domains, Ω_m, all diffusion coefficients are set to zero. In the K⁺ channels, ${\Omega }_{c,{\text{K}}^{+}}$ , the diffusion coefficient for K⁺ is set to $d_{{\text{K}}^{+}}\times D_{{\text{K}}^{+}}$ , where $d_{{\text{K}}^{+}}$ is a channel scaling factor. The diffusion coefficient for the remaining ions are set to zero in the K⁺ channels. Similarly, in the Na⁺ channels, ${\Omega }_{c,{\text{Na}}^{+}}$ , all diffusion coefficients are set to zero when the channel is closed and the diffusion coefficient for Na⁺ is set to the value $d_{{\text{Na}}^{+}}\times D_{{\text{Na}}^{+}}$ when the channel is open. The justification for the choice of scaling factors $d_{{\text{K}}^{+}}$ and $d_{{\text{Na}}^{+}}$ is specified below.