HDL + LDL exchange reactions not balanced with WBM when using custom diet

[Drew Alessi]

Hi Everyone,

I am using the WBM to study the LDL and HDL exchange reactions included in the model. I am using a custom American Average diet adapted from this paper. The diet was created using the VMH.life Diet Designer tool.

My issue is that I noticed that the exchange reactions are not balanced. The sum of my reactions is not 0. I have tested other diets, and they are all at steady-state. 

Does anyone have any suggestions or workarounds to resolve this issue? Please let me know if you require any other information!

Thank you in advance!

Best,

Drew

[Daniel Fässler]

Hi Drew,

Could you provide more information on how you integrated the diet?
What do you mean by saying that the exchange reactions are not balanced? Do you mean that you tried to optimize an objective function, and this resulted in an error because the steady-state assumption does not hold?

Best,

Daniel

[Drew Alessi]

Hi Daniel,

Thank you very much for the assistance!

I am currently using the 'optimizeCBModel' function to solve for all reactions involving 'LDL' and 'HDL' exchange from the [e] to [bc] compartments. I integrated the diet by using the 'SetDietConstraints' function and inputting it into my MATLAB 2019b workspace. After I optimized the model for the whole body objective reaction and recorded the fluxes for each of the reactions involved in HDL-C and LDL-C exchange into [bc], I noticed the sum of the fluxes was not equal to 0. This only occurs with the female WBM. I am using the CPLEX 12.10 solver. Also, I just want to note that I did not receive any infeasibility errors or any other type of error message. I only realized there was an error when I took the sum of the reactions and realized it wasn't 0.

Please let me know if you would like to see my data (list of reactions and their respective fluxes) for additional context!

Best,

Drew

[Daniel Fässler]

Hi Drew,

The fact that you receive no infeasibility error with that function means that the steady-state assumption is fulfilled. There might be a misunderstanding regarding the sum of fluxes: The rate of production equals the rate of consumption for each metabolite. So, the sum of all fluxes for reactions that involve a particular metabolite should equal zero, ensuring that the net flux for that metabolite is zero.  This principle applies to all single metabolites in the network. Why should the recorded fluxes for each of the reactions involved in HDL-C and LDL-C exchange sum up to zero? 

Yes, can you tell me the list of fluxes and reactions that you would except to sum up to zero?

Best,

Daniel

[Drew Alessi]

Hi Daniel,

Thank you again for your help with this issue.  I just want to clarify that I am referring to the steady state mass balance for LDL in the blood compartment [bc]. I agree with everything you mentioned in your last response, and that is why I am confused about the sum of all reactions involving LDL to equal a value of 0.11. I listed the LDL reactions (as well as their reaction equations and respective fluxes) below that I am referring to. I also provided the 'female_solution' structure from MATLAB. All of my code is based on this MATLAB live script, the only thing I changed is the diet. I also provided the diet composition that I used. Please let me know if you would like me to provide any other details or code. I greatly appreciate your time and assistance.

Best,

Drew

Reaction Name | Reaction Formula | Flux Vector
Liver_EX_ldl_hs(e)_[bc] Liver_ldl_hs[e]  <=> ldl_hs[bc]  0.021386
Adipocytes_EX_ldl_hs(e)_[bc] Adipocytes_ldl_hs[e]  <=> ldl_hs[bc]  0.020303
Muscle_EX_ldl_hs(e)_[bc] Muscle_ldl_hs[e]  <=> ldl_hs[bc]  0.024564
Kidney_EX_ldl_hs(e)_[bc] Kidney_ldl_hs[e]  <=> ldl_hs[bc]  -2.7648
Kidney_EX_ldl_hs[bcK]_[bc] Kidney_ldl_hs[bcK]  -> ldl_hs[bc]  2.8093

Sum: 0.110753

solution_female =

  struct with fields:

    origStatText: 'Optimal solution found'
       objLinear: 0
    objQuadratic: 0.0557
               f: 0
              f0: NaN
              f1: NaN
              f2: NaN
               v: [83521×1 double]
               y: [58851×1 double]
               w: [83521×1 double]
               s: [166202×1 double]
          solver: 'ibm_cplex'
            stat: 1
        origStat: 1
            time: 7.9560
          vars_v: []
          ctrs_y: [107351×1 double]
          ctrs_s: [107351×1 double]
               x: [83521×1 double]

[Drew Alessi]

[Daniel Fässler]

Hi Drew,

thank you for your response. Did you check if you forgot a reaction with:

id = find(contains(model.mets, 'ldl_hs[bc]'))

model.rxns(find(model.S(id,:)))

It shows the reactions you provided?
Can you also check:
model.S(id,:)
Are all the entries equal to one?

You are solving a quadratic optimization problem, right? Because in your output there is     objQuadratic: 0.0557 .

Bes regards

Daniel

[Drew Alessi]

Hi Daniel,

Thank you again for your help with this issue.

Yes! We are solving a quadratic optimization problem, the objective function is the minimization of the sum of squares. We are using the code from the Matlab live script and recommended CPLEX solver from this paper.

I ran the codes you provided and unfortunately received the same results as prior. I found the same number of reactions for 'ldl_hs[bc] as previously reported in the male and female models, and all entries are equal to 1.

female.rxns(find(female.S(id,:)))

ans =

  5×1 cell array

    {'Liver_EX_ldl_hs(e)_[bc]'     }
    {'Adipocytes_EX_ldl_hs(e)_[bc]'}
    {'Muscle_EX_ldl_hs(e)_[bc]'    }
    {'Kidney_EX_ldl_hs(e)_[bc]'    }
    {'Kidney_EX_ldl_hs[bcK]_[bc]'  }

female.S(id,:)

ans =

         (1,66567)          1
         (1,69331)          1
         (1,78106)          1
         (1,78225)          1
         (1,78226)          1

Best,

Drew

[Drew Alessi]

Hi Daniel,

Thank you again for helping me work through this problem. Do you have any other suggestions regarding troubleshooting this issue? I greatly appreciate your time.

Best,

Drew

[Daniel Fässler]

Hi Drew,

is it possible that you share your code? I tried to reproduce it but I didn't have the problem that the fluxes do not sum up to zero.

Best,
Daniel

[Drew Alessi]

Hi Daniel,

Thank you for your response and continued assistance! I would be happy to share my code. Just to confirm, are your fluxes summed up to zero when using this custom Average American Diet?

Best,

Drew

[Daniel Fässler]

Hi Drew,

with the diet constraints  you provided, I get a non optimal solution :-(.
I use the optimizeWBModel function:

load Harvetta_1_03c

model = setDietConstraints(female, 'AmericanDiet')
param.minNorm = 1e-6;

% Set no linear objective
model.c = zeros(size(model.c));

% Optimization step
solution = optimizeWBModel(model, param);

Best,
Daniel