%========================================================================== % Please copy paste this file to a MATLAB file with extension *.m % Delete the HTML publishing note (inserted by MATLAB) at the end of this % document %========================================================================== %============================================================================== % Supplement to: Huber, Connolly, Dussmann and Prehn % % Calculations of Figure 3-7 of manuscript and supplementary Figure 1-3 % (control parameter xFig) % % % Includes the following Modules % % - Bioenergetic subsystem "sub_energetic" following % - Beard DA. (2005) PLoS Comp Bio. 1(4), e36 (mitochondrial % respiration) % - Korzeniewski B. (1998) Biophys Chem 75(1), 73-80 (ATP/ADP transfer) % % Heinrich Huber, Royal College of Surgeons in Ireland, 2011 % %========================================================================== % CODE USAGE %========================================================================== % % call regulation(figurename) for calculation of particular figure % % figurename = 'Fig3' for Figure3 % figurename = 'Fig4' for Figure4 % figurename = 'Fig5' for Figure5 % figurename = 'Fig6' for Figure6 % figurename = 'Fig7' for Figure7 % figurename = 'Supl1' for Supplement 1 % figurename = 'Supl2' for Supplement 2 % figurename = 'Supl3' for Supplement 3 % %========================================================================== % S A M P L E C A L L S %========================================================================== % % regulation('Fig3') % % Email comments and suggestions to heinhuber@rcsi.ie %========================================================================== function regulation(xFig) global pp_casp3 Qtot0 Ctot0 W_c cytosolic_ATP global NADtot K_ADTP_dyn c_inc xpar lgpar varpar state_fact t_final txpar % % Control parameter for figure calculation % % % The variable Out_par contains the following elements % % Out_par(1) = Input function (?Dehydrogenase flux?) % Out_par(2) = Electron flux through complex I/II % Out_par(3) = Electron flux through complex III % Out_par(4) = Electron flux through complex IV % Out_par(5) = FoF1 ATP synthase flux % Out_par(6) = ATP transferase (inner membrane) % Out_par(7) = unused % Out_par(8) = unused % Out_par(9) = Matrix Hydrogen % Out_par(10) = Total IM cytochrom-c % Out_par(11) = Total ubiquinol % Out_par(12) = free active caspase3 in nM % Out_par(13) = Mitochondrial Membrane potential % Out_par(14) = Membrame Protomotive Force % Out_par(15) = Hydrogen leaks; % Out_par(16) = cytosolic (external medium) ADP/ATP ratio % Out_par(17) = pH difference mito and cytosol % % Time points for calculation % t_prior = -80; % Time point for prior to steady state t_start = 0; % Time point for steady state t_final = 60; % Time point for prior to steady state t_no_time = t_final+10; % Dummy time for neglected effects F = 0.096484; % kJ mol^{-1} mV^{-1} R = 8314e-6; % Universal gas constant (kJ * mol^{-1} * K^{-1}) T = (273.15 + 37); % Temperature (K), 37 degree RT = R*T; RTF = RT/F; cytosolic_ATP = true; % % % Bioenergetic constants % %************************************************************************* % Model parameters - Compare Figure 2 of main text %************************************************************************* % % Metabolites % xpar = zeros(31,1); xpar(1) = 0.000242; % NADH calculated after Evans_77 xpar(2) = 2.70e-3; % Cyt-c after Beard xpar(3) = 1.35e-3; % total IMS Ubiquinol, Q+QH2, molar xpar(4) = 0.00486; % ATP total Adenosine phosphates in cytosol % calculated from Evans_77 %************************************************************************* % FIG. 2 FLUX 1: Parameter for input function %************************************************************************* xpar(5) = 0.13413e-3; % Dehydrogenase flux input xpar(6) = 0.67668e-3; % Dehydrogenase flux input xpar(7) = 0.09183; % Dehydrogenase activity xpar(8) = 45.807; % Input-flux: Initial disturbance of equilibrium %************************************************************************* % FIG. 2 FLUX 2: Parameter for complex I %************************************************************************* xpar(9) = 1.0200e+003; % Complex I activity %************************************************************************* % FIG. 2 FLUX 3: Parameter for complex III %************************************************************************* xpar(10) = 0.2241; % Complex III activity xpar(11) = 0.19172e-3; % Complex III / Pi parameter xpar(12) = 25.31e-3; % Complex III / Pi parameter %************************************************************************* % FIG. 2 FLUX 4: Parameter for complex IV %************************************************************************* xpar(13) = 3.2131e-004; % Complex IV activity xpar(14) = 1.2e-4; % kinetic constant for complex IV (M) %************************************************************************* % FIG. 2 FLUX 5: Parameter for ATP synthase and Mg-binding to ATP/ADP %************************************************************************* xpar(15) = 6.8294e+003; % F1F0 ATPase activity % % Parameter for cytosolic ATP production and consumption xpar(16) = 0.1980; % kinetic constant, Huber 2011, Materials and Methods xpar(17) = 3; % equilibrium value, Huber 2011, xpar(18) = 0.02; % cytosolic Mg concentration %************************************************************************* % FIG. 2 FLUX 6: Parameter for Adenosine Transferase %************************************************************************* xpar(19) = 0.0020; % ANT activity xpar(20) = 3.50e-6; % ANT Michaelis-Menten constant %************************************************************************* % FIG. 2 FLUX 7: Proton leak activity %************************************************************************* xpar(21) = 150.0; %************************************************************************* % FIG. 2 FLUX 8: Parameter for OM transporters %************************************************************************* xpar(22) = 2e1; % MOM permeability to protons (micron s^{-1}) xpar(23) = 5.99; % MOM area per unit volume micron^{-1} xpar(24) = 85.0; % MOM permeability to phosphate xpar(25) = 327; % MOM permeability to nucleotides %************************************************************************* % FIG. 2 FLUX 9: Phosphate-Hydrogen-Cotransport %************************************************************************* xpar(26) = 10^(-6.75); % H/Pi co-transpor Molar (form factor 1) xpar(27) = 0.45082e-3; % H+/Pi co-transport activity (form factor 2) xpar(28) = 3.3943e5; % H+/Pi? co-transport activity %************************************************************************* % FIG. 2 FLUX 10: Potassium-Hydrogen-Antiport %************************************************************************* xpar(29) = 2.9802e7; % K+ / H+ antiporter activity xpar(30) = 100; % Inner Matrix hydrogen buffer capacity xpar(31) = 6.7568e-6; % Inner Membrane capacitance npars = 31; % Numbers of parameters for sensitivity analysis and PCA % % Symbol definition % c_inc = 1; % Cyt-c after release set to 0.1 % NADtot = xpar(1); % calculated after Evans_77 Ctot0 = xpar(2); % total IMS Cyt-C, Cred+Cox, molar Qtot0 = xpar(3); % total IMS Ubiquinol, Q+QH2, molar ADTP_tot = xpar(4); % total Adenosine phosphates in cytosol, calcuulated from Evans_77 state_3_fact = 100/101; % ATP:ADP ratio = 100:1% (99 % ATP) V_mito = 0.05; % 5 % mitochondrial fraction W_c = 1/V_mito; % cytosolic volume % varpar = false; % for default calculations bioenergetic parameter will be kept fixed % (will be varied for sensitivity analysis and PCA) lgpar = false; % for PCA we have to use variations of logarithmic parameters % i.e. parameters will be transformed to log and then back % % ODE options % txpar = str2mat('Total NADH', 'Total cyt-c','Total ubiquinone','Total adenosine phos.','Input 1','Input 2','Input 3',... 'Input 4','Cmplx I activity','Cmplx III activity','Cmplx III param. 1','Cmplx III param. 2','Cmplx IV activity','Cmplx IV param. 1',... 'ATP Synthase act.','ATP cyt. kin','ATP cyt. equilibrium','Mg concentration ','ANT activity','ANT par. 1','Proton leak activity',... 'MOM proton perm.','MOM area','MOM nucleotide perm.','MOM phosphate perm.','H+/Phos par. 1','H+/Phos par. 2','H+/Phos par. activity',... 'K/H activity','Matrix H+ buffer','Mito capacitance'); % % Calculated selected figure % switch xFig case {'Fig3','Fig4', 'Supl1', 'Supl2', 'Supl3'} % 'Figure 3,4 & Supplement: Single parameter variations' % - Vary each parameter by a factor 0.1, 0.2, 5, 10 % - Calculate changes of model results by parameter variations % - mitochondrial membrane potential in state_3 and state_4 % - pH in state_3 and state_4 % - Respiratory control ratio % - ATP synthase direction and strength after full cyt_c release (0.1% remnant cyt-c) % - ATP synthase direction and strength after partial % cyt_c release (5% remnant cyt-c) %=================================================================================================================================== calc_Fig_3_4 %=================================================================================================================================== case {'Fig5','Fig6','Fig7'} %================================================================== %================================================================= % 'Figure 5, 6 & 7: Parameter Cluster Analysis' % % - Simulate at different cyt-c levels to mimic intact cyt-c or its release % - Analyse change of state variables as consequence of pairwise parameter changes % - Calculate the Hessian matrix (second order approximation) of the error function %=================================================================================================================================== calc_Fig_5_6_7 %=================================================================================================================================== otherwise 'No valid Figure number. ' 'Call: regulation(figurename) ' %=================================================================================================================================== end %======================================================================================== % % Figure subroutines % %======================================================================================== % function calc_Fig_3_4 %=================================================================================================================================== % % - Vary each parameter by a factor 0.1, 0.2, 5, 10 % - Calculate changes of model results by parameter variations % - mitochondrial membrane potential in state_3 and state_4 % - Respiratory control ratio % - ATP synthase direction and strength after full cyt_c release (0.1% remnant cyt-c) % - ATP synthase direction and strength after partial cyt_c release (5% remnant cyt-c) % - Effect of complex I cleavage to mitochondrial membrane potential % - Increase of mitochondrial membrane potential after increase of ATP (0.1% remnant cyt-c) % - Increase of mitochondrial membrane potential after increase of ATP (5% remnant cyt-c) %=================================================================================================================================== t_final = 60; % Time point for prior to steady state pp_casp3 = 0; % Dummy in that case (no active caspase) Ctot0 = 2.70e-3; % Remains unchanged except for parameter variation (no complex I cleavage) Qtot0 = 1.35e-3; % Remains unchanged except for parameter variation (no cytc-c release) K_ADTP_dyn = 0; % Dummy in that case (no cytosolic ATP production) c_inc = 0; % Dummy in that case (no cytc-c release) varpar = true; %Parameters will be varied for sensitivity analysis testmesh = [0.1 0.2 1 5 10]; [idummy,ntestmesh] = size(testmesh); % % Initialise variables % psi_3 = zeros(ntestmesh, npars); % Mitochondrial membrane potential in state 3 (no cytc release) psi_4 = zeros(ntestmesh, npars); % Mitochondrial membrane potential in state 4 (no cytc release) x_control = zeros(ntestmesh, npars); % Respiratory control ratio (no cytc release) psi_feed = zeros(ntestmesh, npars); % Effect of complex I inhibition after cytc release and caspase actvation % f0f1_rel_all = zeros(ntestmesh, npars); % ATP Synthase after full cytc release (0.1 % remnant) f0f1_rel_func = zeros(ntestmesh, npars); % ATP Synthase after partial cytc release (5 % remnant) % psi_rel_all = zeros(ntestmesh, npars); % Mitochondrial membrane potential after full cytc release (0.1 % remnant) psi_rel_func = zeros(ntestmesh, npars); % Mitochondrial membrane potential after partial cytc release (5 % remnant) % % Calculate first steady state as universal starting points for % state-3 and state-4 % Ctot0 = 2.70e-3; Qtot0 = 1.35e-3; state_fact = 3/4; xpar(16) = 1.1; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 ADTP_tot = xpar(4); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0_3 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); %calculate the initial steady state % state_fact = 100/101; % ATP:ADP =100:1 xpar(16) = 0.0540; % kinetic constant, Huber 2011, Materials and Methods (=log(2)/ 3.5 min) xpar(17) = 30; % equilibrium value, Huber 2011, ADTP_tot = xpar(4); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0_4 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); %calculate the initial steady state % Conditions for state-3 state_fact = 3/4; xpar(16) = 1.1; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 state_3 = get_conc_and_fluxes(x0_3); % Conditions for state-4 state_fact = 100/101; % ATP:ADP =100:1 xpar(16) = 0.0540; % kinetic constant, Huber 2011, Materials and Methods (=log(2)/ 3.5 min) xpar(17) = 30; % equilibrium value, Huber 2011, state_4 = get_conc_and_fluxes(x0_4); for igut= 1:npars igut kref = xpar(igut); for jmesh = 1:ntestmesh % % Parameter variation % xpar(igut) = kref*testmesh(jmesh); % ref_c = xpar(2); % Storage for cyt-c ref_q = xpar(3); % Storage for ubiquinone % %=================================================================================================================================== % Test changes of result due to parameter variation %=================================================================================================================================== % Mito membrane potential and control ratio in state-3/4 % % Conditions for state-3 state_fact = 3/4; xpar(16) = 1.1; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 state_3 = get_conc_and_fluxes(x0_3); % Conditions for state-4 state_fact = 100/101; % ATP:ADP =100:1 xpar(16) = 0.0540; % kinetic constant, Huber 2011, Materials and Methods (=log(2)/ 3.5 min) xpar(17) = 30; % equilibrium value, Huber 2011, % Modelled such that glycolysis + mitochondria render ATP:ADP if igut ==16 xpar(16) = 0.0540*testmesh(jmesh); end if igut ==17 xpar(17) = 30*testmesh(jmesh); end state_4 = get_conc_and_fluxes(x0_4); psi_3(jmesh, igut) = state_3(13); % Mitochondrial membrane potential in state_3 psi_4(jmesh, igut) = state_4(13); % Mitochondrial membrane potential in state_4 resp_3(jmesh, igut) = state_3(1); % Respiration in state_3 resp_4(jmesh, igut) = state_4(1); % Respiration in state_4 pmf_3(jmesh, igut) = state_3(12); % Protomotive Force in state_3 pmf_4(jmesh, igut) = state_4(12); % Protomotive Force in state_4 x_control(jmesh, igut) = state_3(1)/state_4(1); f0f1_3(jmesh, igut) = state_3(5); % Protomotive Force in state_3 f0f1_4(jmesh, igut) = state_4(5); % Protomotive Force in state_4 % Test for complete cyt-c release xpar(2) = ref_c /200; % Conditions for state-3 state_fact = 3/4; xpar(16) = 1.1; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 if igut ==16 xpar(16) = 1.1*testmesh(jmesh); end if igut ==17 xpar(17) = 2*testmesh(jmesh); end state_rel_all = get_conc_and_fluxes(x0_3); xpar(2) = ref_c; resp_rel_all(jmesh, igut) = state_rel_all(1); f0f1_rel_all(jmesh, igut) = state_rel_all(5); psi_rel_all(jmesh, igut) = state_rel_all(13); pmf_rel_all(jmesh, igut) = state_rel_all(12); % Test for impaired cyt-c release xpar(2) = ref_c/20; % 5% remnant cyt-c % Conditions for state-3 state_fact = 3/4; xpar(16) = 1.1; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 if igut ==16 xpar(16) = 1.1*testmesh(jmesh); end if igut ==17 xpar(17) = 2*testmesh(jmesh); end state_rel_func = get_conc_and_fluxes(x0_3); xpar(2) = ref_c; resp_rel_func(jmesh, igut) = state_rel_func(1); f0f1_rel_func(jmesh, igut) = state_rel_func(5); psi_rel_func(jmesh, igut) = state_rel_func(13); pmf_rel_func(jmesh, igut) = state_rel_func(12); % Test for complex I cleavage 1% (legacy, not used here) xpar(2) = ref_c /200; xpar(3) = ref_q /100; % Conditions for state-3 state_fact = 3/4; xpar(16) = 1.1; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 if igut ==16 xpar(16) = 1.1*testmesh(jmesh); end if igut ==17 xpar(17) = 2*testmesh(jmesh); end state_rel_feed = get_conc_and_fluxes(x0_3); xpar(2) = ref_c; xpar(3) = ref_q; psi_feed(jmesh, igut) = state_rel_all(13)-state_rel_feed(13); %=================================================================================================================================== % End of test calculation %=================================================================================================================================== end xpar(igut) = kref; end %=================================================================================================================================== % Bar graph output %=================================================================================================================================== ntestp = 20; for i=1:npars test(i) = abs(psi_4(1,i) - psi_4(end,i)); end [dummy, i_psi_4_down ] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp psi_4_reg(i,j) = psi_4(j,i_psi_4_down(i)); end end % for i=1:npars test(i) = abs(psi_3(1,i) - psi_3(end,i)); end [dummy, i_psi_3_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp psi_3_reg(i,j) = psi_3(j,i_psi_3_down(i)); end end for i=1:npars test(i) = abs(resp_3(1,i) - resp_3(end,i)); end [dummy, i_resp_3_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp resp_3_reg(i,j) = resp_3(j,i_resp_3_down(i)); end end % for i=1:npars test(i) = abs(resp_4(1,i) - resp_4(end,i)); end [dummy, i_resp_4_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp resp_4_reg(i,j) = resp_4(j,i_resp_4_down(i)); end end for i=1:npars test(i) = abs(pmf_3(1,i) - pmf_3(end,i)); end [dummy, i_pmf_3_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp pmf_3_reg(i,j) = pmf_3(j,i_pmf_3_down(i)); end end % for i=1:npars test(i) = abs(pmf_4(1,i) - pmf_4(end,i)); end [dummy, i_pmf_4_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp pmf_4_reg(i,j) = pmf_4(j,i_pmf_4_down(i)); end end for i=1:npars test(i) = abs(f0f1_3(1,i) - f0f1_3(end,i)); end [dummy, i_f0f1_3_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp f0f1_3_reg(i,j) = f0f1_3(j,i_f0f1_3_down(i)); end end % for i=1:npars test(i) = abs(f0f1_4(1,i) - f0f1_4(end,i)); end [dummy, i_f0f1_4_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp f0f1_4_reg(i,j) = f0f1_4(j,i_f0f1_4_down(i)); end end for i=1:npars test(i) = abs(x_control(1,i) - x_control(end,i)); end [dummy, i_x_control_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp x_control_reg(i,j) = x_control(j,i_x_control_down(i)); end end % for i=1:npars test(i) = abs(psi_rel_func(1,i) - psi_rel_func(end,i)); end [dummy, i_psi_rel_func_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp psi_rel_func_reg(i,j) = psi_rel_func(j,i_psi_rel_func_down(i)); end end % for i=1:npars test(i) = abs(psi_rel_all(1,i) - psi_rel_all(end,i)); end [dummy, i_psi_rel_all_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp psi_rel_all_reg(i,j) = psi_rel_all(j,i_psi_rel_all_down(i)); end end % % for i=1:npars test(i) = abs(resp_rel_func(1,i) - resp_rel_func(end,i)); end [dummy, i_resp_rel_func_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp resp_rel_func_reg(i,j) = resp_rel_func(j,i_resp_rel_func_down(i)); end end % for i=1:npars test(i) = abs(resp_rel_all(1,i) - resp_rel_all(end,i)); end [dummy, i_resp_rel_all_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp resp_rel_all_reg(i,j) = resp_rel_all(j,i_resp_rel_all_down(i)); end end % % for i=1:npars test(i) = abs(pmf_rel_func(1,i) - pmf_rel_func(end,i)); end [dummy, i_pmf_rel_func_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp pmf_rel_func_reg(i,j) = pmf_rel_func(j,i_pmf_rel_func_down(i)); end end % for i=1:npars test(i) = abs(pmf_rel_all(1,i) - pmf_rel_all(end,i)); end [dummy, i_pmf_rel_all_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp pmf_rel_all_reg(i,j) = pmf_rel_all(j,i_pmf_rel_all_down(i)); end end % % for i=1:npars test(i) = abs(f0f1_rel_func(1,i) - f0f1_rel_func(end,i)); end [dummy, i_f0f1_rel_func_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp f0f1_rel_func_reg(i,j) = f0f1_rel_func(j,i_f0f1_rel_func_down(i)); end end % for i=1:npars test(i) = abs(f0f1_rel_all(1,i) - f0f1_rel_all(end,i)); end [dummy, i_f0f1_rel_all_down] = sort(test,2,'descend'); for j = 1:ntestmesh for i = 1:ntestp f0f1_rel_all_reg(i,j) = f0f1_rel_all(j,i_f0f1_rel_all_down(i)); end end % % Figure 3: State-3 / state-4 mito potential % if strcmp (xFig, 'Fig3') figure subplot(2,2,1), bar(psi_4_reg(1:ntestp,:),'group') axis([0 ntestp+1 100 200]) title('Figure 3A: Mitochondrial Membrane Potential in state-4') subplot(2,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_psi_4_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(2,2,3), bar(psi_3_reg(1:ntestp,:),'group'); axis([0 ntestp+1 100 180]) title('Figure 3B: Mitochondrial Membrane Potential in state-3') subplot(2,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_psi_3_down(1:ntestp),:)),'HorizontalAlignment','left') end if strcmp (xFig, 'Supl1') % % Respiration in state-3 and state-4 % figure subplot(3,2,1), bar(resp_4_reg(1:ntestp,:),'group') title('Supplementary Figure 1A: Respiration in state-4') subplot(3,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_resp_4_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(3,2,3), bar(resp_3_reg(1:ntestp,:),'group'); title('Supplementary Figure 1B: Respiration in state-3') subplot(3,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_resp_3_down(1:ntestp),:)),'HorizontalAlignment','left') end % % Protomotive Force in state-3 and state-4 % if strcmp (xFig, 'Supl2') figure subplot(3,2,1), bar(pmf_4_reg(1:ntestp,:),'group') title('Supplementary Figure 2A: Protomotive Force in state-4') subplot(3,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_pmf_4_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(3,2,3), bar(pmf_3_reg(1:ntestp,:),'group'); title('Supplementary Figure 2B: Protomotive Force in state-3') subplot(3,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_pmf_3_down(1:ntestp),:)),'HorizontalAlignment','left') end % % ATP synthase in state-3 and state-4 % if strcmp (xFig, 'Supl3') figure subplot(3,2,1), bar(f0f1_4_reg(1:ntestp,:),'group') title('Supplementary Figure 3A: Protomotive Force in state-4') subplot(3,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_f0f1_4_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(3,2,3), bar(f0f1_3_reg(1:ntestp,:),'group'); title('Supplementary Figure 3B: Protomotive Force in state-3') subplot(3,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_f0f1_3_down(1:ntestp),:)),'HorizontalAlignment','left') end % % Depolarisation after cyt-c release % if strcmp (xFig, 'Fig4') figure subplot(2,2,1), bar(psi_rel_all_reg(1:ntestp,:),'group'); axis([0 ntestp+1 20 150]) title('Figure 4A: Mito depolarisation after release to 0.1% cyt-c and intact complex I/II') subplot(2,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_psi_rel_all_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(2,2,3), bar(psi_rel_func_reg(1:ntestp,:),'group'); axis([0 ntestp+1 20 150]) title('Figure 4B: Mito depolarisation after release to 5% cyt-c and intact complex I/II') subplot(2,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_psi_rel_func_down(1:ntestp),:)),'HorizontalAlignment','left') end % % % Respiration after cyt-c release % if strcmp (xFig, 'Supl1') figure subplot(2,2,1), bar(resp_rel_all_reg(1:ntestp,:),'group'); title('Supplementary Figure 1C: Respiration after release to 0.1% cyt-c ') subplot(2,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_resp_rel_all_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(2,2,3), bar(resp_rel_func_reg(1:ntestp,:),'group'); title('Supplementary Figure 1D: Respiration after release to 5% cyt-c ') subplot(2,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_resp_rel_func_down(1:ntestp),:)),'HorizontalAlignment','left') end % % Protomotive Force after cyt-c release % if strcmp (xFig, 'Supl2') figure subplot(2,2,1), bar(pmf_rel_all_reg(1:ntestp,:),'group'); title('Supplementary Figure 2C: Protomotive Force after release to 0.1% cyt-c ') subplot(2,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_pmf_rel_all_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(2,2,3), bar(pmf_rel_func_reg(1:ntestp,:),'group'); title('Supplementary Figure 2D: Protomotive Force after release to 5% cyt-c ') subplot(2,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_pmf_rel_func_down(1:ntestp),:)),'HorizontalAlignment','left') end % % ATP synthase after cyt-c release % if strcmp (xFig, 'Supl3') figure subplot(2,2,1), bar(f0f1_rel_all_reg(1:ntestp,:),'group'); axis([0 ntestp+1 -3.e-4 2e-4]) title('Supplementary Figure 4C: ATP synthase after release to 0.1% cyt-c') subplot(2,2,2), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_f0f1_rel_all_down(1:ntestp),:)),'HorizontalAlignment','left') subplot(2,2,3), bar(f0f1_rel_func_reg(1:ntestp,:),'group'); axis([0 ntestp+1 -1e-3 2e-3]) title('Supplementary Figure 4D: ATP synthase after release to 5% cyt-c') subplot(2,2,4), bar(0,'group', 'FaceColor', 'b'); hold on axis([0 2 0 2]) text(0.1,1,strcat(txpar(i_f0f1_rel_func_down(1:ntestp),:)),'HorizontalAlignment','left') end end % % function calc_Fig_5_6_7 %=================================================================================================================================== % % 'Figure 5, 6 and 7 Sloppy parameter analysis according Gutentkunst with PCA' % % - Analyse change of state variables as consequence of pairwise parameter changes % - Calculate the Hessian matrix (second order approximation) % of the error function (i.e. change of state varioiable according parameter variation) %=================================================================================================================================== t_final = 60; % Time point for prior to steady state pp_casp3 = 0; % Dummy in that case (no active caspase) Qtot0 = 0; % Dummy in that case (no complex I cleavage) Ctot0 = 0; % Dummy in that case (no cytc-c release) K_ADTP_dyn = 0; % Dummy in that case (no cytosolic ATP production) c_inc = 0; % Dummy in that case (no cytc-c release) varpar = true; %Parameters will be varied for sloppyness analysis n_state = 1; n_cyt = 1; switch xFig case 'Fig5' deg_cyct(1) = 1; % Conditions for state-3 state_fact = 0.74; % note 0.74 instead of 0.75 was used due to numeric instabilities state_3_fact(1) = 0.74; xpar(16) = 1.2; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 case 'Fig6' state_3_fact(1) = 100/101; deg_cyct(1) = 1; % Conditions for state-4 state_fact = 100/101; % ATP:ADP =100:1 xpar(16) = 0.0540; % kinetic constant, Huber 2011, Materials and Methods (=log(2)/ 3.5 min) xpar(17) = 30; % equilibrium value, Huber 2011, case 'Fig7' state_3_fact(1) = 0.74; deg_cyct(1) = 0.05; % Conditions for state-3 state_fact = 0.74; % note 0.74 instead of 0.75 was used due to numeric instabilities xpar(16) = 1.2; % kinetic constant, (=log(2)/ 0.5 min) xpar(17) = 2; % Modelled so that glycolysis + mitochondria render ATP:ADP = 3:1 end npars = 31; % Numbers of varied parameters xpar = log(xpar); lgpar = true; dpar = 0.002; % best parameter over rage E-2-E-4 If too large:not accurate, too low: numeric noise) % % Legacy: The following two loops (n,i) are only calculated once % for n = 1:n_cyt xpar(2) = log(2.70e-3*deg_cyct(n)); % total IMS Cyt-C, Cred+Cox, molar for k = 1:n_state % sum over all energetic conditions Hessian = zeros(npars,npars); state_fact = state_3_fact(k); % get actual ATP ADTP_tot = exp(xpar(4)); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); %calculate the initial steady state x0mem = x0; % for i = 1:npars % sum over all Gutenkunst variables in pairs for j = 1:i xmem1 = xpar(i); % store variable xmem2 = xpar(j); % store variable if i == 4 % If total ATP levels are changed the initial state has to be recalculated % The subsequent calculations getconc only shifts ATP and ADP, but not total levels % xpar(i) = xmem1 + dpar; % infinitesimal change of parameter for centre derivative ADTP_tot = exp(xpar(i)); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); % calculate the initial steady state Sub_p1 = get_conc(x0); xpar(i) = xmem1 - dpar; % infinitesimal change of parameter for centre derivative ADTP_tot = exp(xpar(i)); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); %calculate the initial steady state Sub_m1 = get_conc(x0); xpar(i) = xmem1; % retrieve variable x0 = x0mem; % xpar(j) = xmem2 + dpar; % infinitesimal change of parameter for centre derivative Sub_p2 = get_conc(x0); xpar(j) = xmem2 - dpar; % infinitesimal change of parameter for centre derivative Sub_m2 = get_conc(x0); elseif j==4 xpar(i) = xmem1 + dpar; % infinitesimal change of parameter for centre derivative Sub_p1 = get_conc(x0); xpar(i) = xmem1 - dpar; % infinitesimal change of parameter for centre derivative Sub_m1 = get_conc(x0); xpar(i) = xmem1; % retrieve variable xpar(j) = xmem2 + dpar; % infinitesimal change of parameter for centre derivative ADTP_tot = exp(xpar(j)); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); % calculate the initial steady state Sub_p2 = get_conc(x0); xpar(j) = xmem2 - dpar; % infinitesimal change of parameter for centre derivative ADTP_tot = exp(xpar(j)); % total Adenosine phosphates in cytosol, calculated from Evans_77 ATP_e = state_fact*ADTP_tot; % ATP level ADP_e = ADTP_tot-ATP_e; % ADP level x0 = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0); %calculate the initial steady state Sub_m2 = get_conc(x0); x0 = x0mem; else xpar(i) = xmem1 + dpar; % infinitesimal change of parameter for centre derivative Sub_p1 = get_conc(x0); xpar(i) = xmem1 - dpar; % infinitesimal change of parameter for centre derivative Sub_m1 = get_conc(x0); xpar(i) = xmem1; % retrieve variable xpar(j) = xmem2 + dpar; % infinitesimal change of parameter for centre derivative Sub_p2 = get_conc(x0); xpar(j) = xmem2 - dpar; % infinitesimal change of parameter for centre derivative Sub_m2 = get_conc(x0); end for m = 1:20 if Sub_p1(m) ~=0 diff1(m) = (Sub_p1(m) - Sub_m1(m))/Sub_p1(m); % Calculate derivative (logarithmic parameter) else diff1(m) = 0; end if Sub_p1(m) ~=0 diff2(m) = (Sub_p2(m) - Sub_m2(m))/Sub_p2(m); % Calculate derivative (logarithmic parameter) else diff2(m) = 0; end end Hessian(i,j) = Hessian(i,j) +diff1*diff2'/(dpar*dpar*4); xpar(i) = xmem1; % set back to original value xpar(j) = xmem2; % set back to original value x0 = x0mem; % back to defauls initial state without change of ATP levels end end for i = 1:npars for j = 1:i Hessian(j,i) = Hessian(i,j); % Hessian is symmetric end end %************************************************************************* % Calculate Eigenvectors %************************************************************************* [coeff2, latent2, expl2] = pcacov(Hessian); coeff2(:,npars+1) = coeff2(:,npars); coeff2(npars+1,:) = coeff2(npars,:); figure surf(-coeff2(:,1:5));view(0,-90) xlabel('4 major Principal Components') ylabel('Original Parameters, Suppl. Table VII') zlabel('Parameter Loading') set (gca, 'LineWidth',2) 'Explanatory factors of principal components' expl2 latentmem(n+(k-1)*2,:,:) = coeff2(:,:); expl2mem(n+(k-1)*2,:) = expl2(:); var(n+(k-1)*2,1) = state_3_fact(n); var(n+(k-1)*2,2) = state_3_fact(k); if cytosolic_ATP ipargroup = [1 2 3 9 10 11 12 13 14 0 0 0 0 % Complex Metabolite and Parameters 4 15 18 19 20 26 27 28 0 0 0 0 0 % ATP production 16 17 0 0 0 0 0 0 0 0 0 0 0 % cytosolic ATP 5 6 7 8 21 22 23 24 25 29 30 31 0]; % Protein leaks & Rest np(1) = 9; np(2) = 8; np(3) = 2; np(4) = 12; figure for i = 1:3 for j = 1:4 sum = 0; for k = 1:np(j) sum = sum + coeff2(ipargroup(j,k),i)*coeff2(ipargroup(j,k),i); end pc(i,j) = sum; end subplot(2,2,i),pie(pc(i,:)) title (strcat('Cluster contributions, Cluster: ',num2str(i))) end subplot(2,2,4), pie([1/10 2/10 3/10 4/10 ]) title('Dummy Figure for legend') legend('Complex Metabolite & Parameters','ATP production','cytosolic ATP','Protein leaks & Rest') end end end % save ('latentmem') end % % end function [Init]=initial(ADP_e, ATP_e, Ctot0, Qtot0) % % Sets the initial conditions % Init(1) = 10^(-8); % Matrix hydrogen concentration pH=8 H_x Init(2) = 0.050; % Matrix potassium concentration [50 mM] K_x Init(3) = 0.001; % Matrix magnesium concentration [1 mM] Mg_x Init(4) = 1.5e-3; % Reduced NADH in matrix [1.5mM] NADH_x Init(5) = 0.0008; % Reduced ubiquinol in matrix [0.8 mM] QH2 Init(6) = 1.0e-3; % Reduced cyt-c in IMS [1 mM] Cred Init(7) = 2.6e-5; % O2 unsued as state variable, parameter only Init(8) = ATP_e; % Matrix free ATP (set to external as initial) ATP_x Init(9) = ADP_e; % Matrix free ADP (set to external as initial) ADP_x Init(10) = 0; % Matrix ATP bound to magnesium ATP_mx Init(11) = 0; % Matrix ADP bound to magnesium ADP_mx Init(12) = 10e-3; % Matrix inorganic phosphate [ 10 mM] Pi_x Init(13) = 0; % IMS free ATP ATP_i Init(14) = 0; % IMS free ADP ADP_i Init(15) = 0; % IMS free AMP AMP_i Init(16) = 0; % IMS ATP bound to magnesium Init(17) = 0; % IMS ADP bound to magnesium Init(18) = 10e-3; % Matrix inorganic phosphate [ 10 mM] Pi_i Init(19) = 150; % Mitochondrial membrane potential [150 mV] dPsi Init(20) = Ctot0; % total cyt-c in IMS Ctot0 Init(21) = Qtot0; % total ubiquinol Qtot0 Init(22) = 10^(-8.0); % IMS hydrogen concentration pH=8 H_i Init(23) = ATP_e; % cytosolic ATP (given according text) ATP_c Init(24) = ADP_e; % cytosolic ADP (given according text) ADP_c end % % Subroutines for Sensitivity and PCA % function [initial_equilibrium] = get_first_steady_state(ADP_e, ATP_e, Ctot0, Qtot0) % % ODE options % options = odeset('RelTol',1e-5, 'AbsTol',1e-8, 'MaxStep',10e-1, ... 'InitialStep',1e-1, 'MaxOrder',5, 'BDF','on'); % t_prior = -80; % Time point for prior to steady state t_start = 0; % Time point for steady state % xo = initial(ADP_e, ATP_e, Ctot0, Qtot0); % calculate initial conditions % % calculate steady state % [t0_vec, x_steady]= ode15s(@sub_energetic,[t_prior t_start],xo,options); % % calculate time series of state transition for isolated mitochondria % initial_equilibrium = x_steady(end,:)'; end function [final_equilibrium] = get_conc(initial_equilibrium); global t_final % % ODE options % options = odeset('RelTol',1e-5, 'AbsTol',1e-8, 'MaxStep',10e-1, ... 'InitialStep',1e-1, 'MaxOrder',5, 'BDF','on'); % t_start = 0; % Time point for steady state % [t_out,x1_mat] = ode15s(@sub_energetic,[t_start t_final],initial_equilibrium,options); % % return vector of steady state % final_equilibrium = x1_mat(end,:)'; end function [final_equilibrium] = get_conc_and_fluxes(initial_equilibrium); global C t_final % % ODE options % options = odeset('RelTol',1e-5, 'AbsTol',1e-8, 'MaxStep',10e-1, ... 'InitialStep',1e-1, 'MaxOrder',5, 'BDF','on'); % t_start = 0; % Time point for steady state % [t_out,x1_mat] = ode15s(@sub_energetic,[t_start t_final],initial_equilibrium,options); % % states_and_fluxes = []; for i = 1:length(t_out); t = t_out(i); x = x1_mat(i,:)'; f = sub_energetic(t,x); states_and_fluxes = [states_and_fluxes;C]; end states_and_fluxes(1,end) = states_and_fluxes(1,end)/states_and_fluxes(1,1); % normalize 1st (respiration flux) component to 100% initial % % return vector of steady state % final_equilibrium = states_and_fluxes(end,:)'; end %======================================================================================== % % E N E R G E T I C S U B - R O U T I N E % %======================================================================================== function [f] = sub_energetic(t,x); %======================================================================================== % Code supplement to Huber, Duessmann, Rehm and Prehn % % Bioenergetic part of the model including cytosolic extensions % Heinrich Huber, Royal College of Surgeons in Ireland, 2010 % % Email comments and suggestions to heinhuber@rcsi.ie %========================================================================== %========================================================================== %============== % The bioenergetic model is derived from the great work of % % Daniel A. Beard, Department of Physiology, Medical College of Wisonsin, % Milwaukee, WI. % % Beard DA. (2005) A Biophysical Model of the Mitochondrial Respiratory % System and Oxidative Phosphorylation. PLoS Comp Bio. 1 (4) e36 % % Including the ATP/ADP transfer modelling from: % Korzeniewski B. (1998), Biophys Chem 75(1), 73-80 %======================================================================================== % % We have included the following features % % - Cyt-c release for the total and for one redox fraction % % - Facultative feedback of Caspase-3 to Complex I/II % % - Facultative blockage of ATP synthase % % - Cytosolic ATP production and consumption % %========================================================================== %============== %======================================================================================== % % Notes % % - For consistency and comparability to Beard's original model % we kept the constants in seconds and then rescaled them to minutes % % - To maintain redox equilibrium, always the total and one redox % fraction of cyt-c/ubiqinon is affected from release/cleavage % % - Complex II is not explicitetly considered (similar to Beard). % Caspase-3 has been reported to cleave both, complex I and II. % Likewise, as they work in parallel, they are lumped together %======================================================================================== % % % INPUTS: % t - time % x - vector of state variables % % OUTPUTS: % f - vector of time derivatives of concentrations % % Global INPUT % pp_casp3: Spline parameter of temporal concentrations for free, active % caspase-3 (in uM) starting from time of depolarization (i.e. % shifted) % W_c: Cytosolic volume fraction % Qtot0: Initial ubiquinon (M) % Ctot0: Initial cyt-c in intermembrane space (M) % % Global OUTPUT % C - Output vector for plotted concentration % %======================================================================================== % global C cytosolic_ATP global NADtot xpar varpar lgpar state_fact ADTP_tot2 W_c % %************************************************************************* % Concentration state variables: %************************************************************************* H_x = x(1); % Matrix hydrogen concentration K_x = x(2); % Matrix potassium concentration Mg_x = x(3); % Matrix magnesium concentration NADH_x = x(4); % Reduced NADH in matrix QH2 = x(5); % Reduced ubiquinol in matrix Cred = x(6); % Reduced cyt-c in IMS O2 = x(7); % O2 unsued as state variable, parameter only ATP_x = x(8); % Matrix free ATP (set to external as initial) ADP_x = x(9); % Matrix free ADP (set to external as initial) ATP_mx = x(10); % Matrix ATP bound to magnesium ADP_mx = x(11); % Matrix ADP bound to magnesium Pi_x = x(12); % Matrix inorganic phosphate ATP_i = x(13); % IMS free ATP ADP_i = x(14); % IMS free ADP AMP_i = x(15); % IMS free AMP ATP_mi = x(16); % IMS ATP bound to magnesium ADP_mi = x(17); % IMS ADP bound to magnesium Pi_i = x(18); % Matrix inorganic phosphate dPsi = x(19); % Mitochondrial membrane potential Ctot = x(20); % total cyt-c in IMS Qtot = x(21); % total ubiquinol H_i = x(22); % IMS hydrogen concentration ATP_e = x(23); % cytosolic ATP ADP_e = x(24); % cytosolic ADP AMP_e = 0; % cytosolic AMP concentration (Molar), Nicholls and others % %************************************************************************* % Model parameter %************************************************************************* % % Thermodynamic parameter % F = 0.096484; % kJ mol^{-1} mV^{-1} R = 8314e-6; % Universal gas constant (kJ * mol^{-1} * K^{-1}) T = (273.15 + 37); % Temperature (K), 37 degree RT = R*T; % % Parameters for literature sources, please refer to Beard PLoS Comp Bio. 1 (4) e36) % dG_C1o = -69.37; % kJ mol^{-1} dG_C3o = -32.53; % kJ mol^{-1} dG_C4o = -122.94; % kJ mol^{-1} dG_F1o = 36.03; % kJ mol^{-1} n_A = 3.0; % numbers of proteins used by ATP synthase % % Cytosolic Ion concentration and pH % pH_e = 7.4; % External pH (cytosol) H_e = 10^(-pH_e); % cytosolic hydrogen concentration (Molar) K_ei = 120.0e-3; % cytosolic and IM potassium-concentration (Nicholls, Bioenergetics 2) Pi_e = 0.02; % cytosolic and IM phosphate-concentration % % % Constant for Mg-binding to ATP/ADP % K_DT = 24e-6; % Mg/ATP binding constant M K_DD = 347e-6; % Mg/ADP binding constant M x_MgA = 1e6; % Mg2+ binding activity % if lgpar xpar1 = exp(xpar); % for PCA logarithmic parameter have been used, needed to be back transformed else xpar1 = xpar; % default case end if varpar % During the sensitivity analysis and PCA, below parameter will be changed in the main programme % Likewise, x(23) and x(24) do not change over time NADtot = xpar1(1); % calculated after Evans_77 Ctot = xpar1(2); % total cyt-c in IMS Qtot = xpar1(3); % total ubiquinol ADTP_tot= xpar1(4); % total Adenosine phosphates in cytosol, calculated from Evans_77 % ATP_e = state_fact*ADTP_tot; % State 3 ATP level % ADP_e = ADTP_tot-ATP_e; % State 3 ADP level else % defauilt case for all figures in the main text end if cytosolic_ATP else ATP_e = state_fact*ADTP_tot; % fixed ATP for isolated mitochondria ADP_e = ADTP_tot-ATP_e; % fixed ATP for isolated mitochondria end % % Parameter for input function % k_Pi1 = xpar1(5); % Dehydrogenase flux input k_Pi2 = xpar1(6); % Dehydrogenase flux input x_DH = xpar1(7); % Dehydrogenase activity r_DH = xpar1(8); % Input-flux: Initial disturbance of equilibrium % % Parameter for complex I % x_C1 = xpar1(9); % Complex I activity % % Parameter for complex III % x_C3 = xpar1(10); % Complex III activity k_Pi3 = xpar1(11); % Complex III / Pi parameter k_Pi4 = xpar1(12); % Complex III / Pi parameter % % Parameter for complex IV % x_C4 = xpar1(13); % Complex IV activity k_O2 = xpar1(14); % kinetic constant for complex IV (M) % % Parameter for ATP synthase % x_F1 = xpar1(15); % F1F0 ATPase activity % % Parameter for cytosolic ATP production and consumption % x_ATPK = xpar1(16); % kinetic constant to establish equilibrium K_ADTP_dyn = xpar1(17); % Equilibrium value Mg_tot = xpar1(18); % cytosolic and IM magnesium-concentration (Nicholls, Bioenergetics 2). % % Parameter for Adenosine Transferase % x_ANT = xpar1(19); % ANT activity k_mADP = xpar1(20); % ANT Michaelis-Menten constant % % Proton leaks % x_Hle = xpar1(21); % Proton leak activity % % Parameter for OM transporters % x_Ht = xpar1(22); % mito outer membrane permeability to protons (micron s^{-1}) gamma = xpar1(23); % mito outer membrane area per unit volume micron^{-1} x_A = xpar1(24); % mito outer membrane permeability to nucleotides x_Pi2 = xpar1(25); % mito outer membrane permeability to phosphate % % Phosphate-Hydrogen-Cotransport % k_dHPi = xpar1(26); % H/Pi binding constant Molar k_PiH = xpar1(27); % H+ / Pi ? co-transport Michaelis constant x_Pi1 = xpar1(28); % H+ / Pi ? co-transport activity % % Potassium-Hydrogen-Antiport % x_KH = xpar1(29); % K+ / H+ antiporter activity % % Matrix buffer and membrane capacitance % x_buff = xpar1(30); % Inner Matrix hydrogen buffer capacity CIM = xpar1(31); % Inner Membrane capacitance % % Release and compartment parameter % W_m = 0.143/0.20; % mitochondrial water space (ml water / ml mito) W_x = 0.9*W_m; % Matrix water space (ml water / ml mito) W_i = 0.1*W_m; % IM water space (ml water / ml mito) % IMS water space with the rest of the cell (0.1%) % c_inc is used to modify OXPHOS-contributing cyt-c after release %************************************************************************* % potassium uniporter and adenylate kinase neglected %************************************************************************* x_K = 0; % Passive potassium transporter activity x_AK = 0e6; % AK activity K_AK = 0; % Adenelyte Kinase switched off %************************************************************************* % Balancing moeities %************************************************************************* % NAD_x = NADtot - NADH_x; Q = Qtot - QH2; Cox = Ctot - Cred; null = 0; % %************************************************************************* % Calculate %************************************************************************* % ATP_fx = ATP_x - ATP_mx; ADP_fx = ADP_x - ADP_mx; ATP_fi = ATP_i - ATP_mi; ADP_fi = ADP_i - ADP_mi; % %************************************************************************* % Preventing numerical instabilities %************************************************************************* % Q = Q*(Q>0); QH2 = QH2*(QH2>0); Qtot = Qtot*(Qtot>0); Cox = Cox*(Cox>0); Cred = Cred*(Cred>0); Ctot = Ctot*(Ctot>0); ATP_fx = ATP_fx*(ATP_fx>0); ATP_fi = ATP_fi*(ATP_fi>0); ADP_fx = ADP_fx*(ADP_fx>0); ADP_fi = ADP_fi*(ADP_fi>0); % %************************************************************************* % ADP/Mg/K binding in E space %************************************************************************* ADP_me = ( (K_DD+ADP_e+Mg_tot) - sqrt((K_DD+ADP_e+Mg_tot)^2-4*(Mg_tot*ADP_e)) )/2; ADP_fe = ADP_e - ADP_me; Mg_e = Mg_tot - ADP_me; Mg_i = Mg_e; % Mg in inner space same as external K_i = K_ei; % % %************************************************************************* % Calculating Membrane proton motive force and respiration fluxes %************************************************************************* dG_H = F*dPsi + 1*RT*log(H_i/H_x); % Protomotive force dG_C1op = dG_C1o - 1*RT*log(H_x/1e-7); dG_C3op = dG_C3o + 2*RT*log(H_x/1e-7); dG_C4op = dG_C4o - 2*RT*log(H_x/1e-7); dG_F1op = dG_F1o - 1*RT*log(H_x/1e-7); % J_DH = x_DH*(r_DH*NAD_x-NADH_x)*((1+Pi_x/k_Pi1)/(1+Pi_x/k_Pi2)); J_C1 = x_C1*(exp(-(dG_C1op+4*dG_H)/RT)*NADH_x*Q - NAD_x*QH2); J_C3 = x_C3*((1+Pi_x/k_Pi3)/(1+Pi_x/k_Pi4))*... (exp(-(dG_C3op+4*dG_H-2*F*dPsi)/(2*RT))*Cox*QH2^0.5 - Cred*Q^0.5); J_C4 = x_C4*(O2/(O2+k_O2))*(Cred/Ctot)*... (exp(-(dG_C4op+2*dG_H)/(2*RT))*Cred*(O2^0.25) - Cox*exp(F*dPsi/RT)); J_F1 = x_F1*(exp(-(dG_F1op-n_A*dG_H)/RT)*(K_DD/K_DT)*ADP_mx*Pi_x - ATP_mx); % %************************************************************************* % Modelling ATP transferase Korzeniewski 1998 %************************************************************************* Psi_x = -0.65*dPsi; Psi_i = +0.35*dPsi; if (ADP_fi > null) || (ATP_fi > null) J_ANT = x_ANT*( ADP_fi/(ADP_fi+ATP_fi*exp(-F*Psi_i/RT)) - ADP_fx/(ADP_fx+ATP_fx*exp(-F*Psi_x/RT)) )*(ADP_fi/(ADP_fi+k_mADP)); else J_ANT = 0; end % %************************************************************************* % Calculating ionic fluxes %************************************************************************* % H2PIi = Pi_i*H_i/(H_i+k_dHPi); H2PIx = Pi_x*H_x/(H_x+k_dHPi); J_Pi1 = x_Pi1*(H_x*H2PIi - H_i*H2PIx)/(H2PIi+k_PiH); J_Hle = x_Hle*dPsi*(H_i*exp(F*dPsi/RT)-H_x)/(exp(F*dPsi/RT)-1); J_KH = x_KH*( K_i*H_x - K_x*H_i ); J_K = x_K*dPsi*(K_i*exp(F*dPsi/RT)-K_x)/(exp(F*dPsi/RT)-1); J_AKi = x_AK*( K_AK*ADP_i*ADP_i - AMP_i*ATP_i ); J_AMP = gamma*x_A*(AMP_e-AMP_i); J_ADP = gamma*x_A*(ADP_e-ADP_i); J_ATP = gamma*x_A*(ATP_e-ATP_i); J_Pi2 = gamma*x_Pi2*(Pi_e-Pi_i); J_Ht = gamma*x_Ht*(H_e-H_i); % J_MgATPx = x_MgA*(ATP_fx*Mg_x-K_DT*ATP_mx); J_MgADPx = x_MgA*(ADP_fx*Mg_x-K_DD*ADP_mx); J_MgATPi = x_MgA*(ATP_fi*Mg_i-K_DT*ATP_mi); J_MgADPi = x_MgA*(ADP_fi*Mg_i-K_DD*ADP_mi); % %************************************************************************* % Cytosolic Energy balance for single cell model (see text) %************************************************************************* J_ATPK = x_ATPK * (ATP_e- K_ADTP_dyn*ADP_e); % %************************************************************************* % Calculating derivatives for next step %************************************************************************* % f(1) = x_buff*H_x*( +1*J_DH - (4+1)*J_C1 - (4-2)*J_C3 - (2+2)*J_C4 + (n_A-1)*J_F1 + 2*J_Pi1 + J_Hle - J_KH )/W_x; % H_x f(2) = (J_KH + J_K)/W_x; % K_x f(3) = (-J_MgATPx - J_MgADPx)/W_x; % Mg_x f(4) = (+J_DH - J_C1)/W_x; % NADH f(5) = (+J_C1 - J_C3)/W_x; % QH2 f(6) = (+2*J_C3 - 2*J_C4)/W_i; % Include loss of reduce cyt-c with half time t12cyto f(7) = 0; % O2 unused f(8) = (+J_F1 - J_ANT)/W_x; % ATP_x f(9) = (-J_F1 + J_ANT)/W_x; % ADP_x f(10) = (J_MgATPx)/W_x; % ATP_mx f(11) = (J_MgADPx)/W_x; % ADP_mx f(12) = (-J_F1 + J_Pi1 )/W_x; % Pi_x MOD f(13) = (+J_ATP + J_ANT + J_AKi )/W_i; % ATP_i f(14) = (+J_ADP - J_ANT - 2*J_AKi )/W_i; % ADP_i f(15) = (+J_AMP + J_AKi)/W_i; % AMP_i f(16) = (J_MgATPi)/W_i; % ATP_mi f(17) = (J_MgADPi)/W_i; % ADP_mi f(18) = (-J_Pi1 + J_Pi2 )/W_i; % Pi_i f(19) = ( 4*J_C1 + 2*J_C3 + 4*J_C4 - n_A*J_F1 - J_ANT - J_Hle - J_K )/CIM; %MOD f(20) = 0; % unused f(21) = 0; % unused f(22) = (- J_DH +(4+1)*J_C1 + (4-2)*J_C3 + (2+2)*J_C4 - (n_A-1)*J_F1 ... - 2*J_Pi1 - J_Hle + J_KH + J_Ht)/W_i; % Change Intermembrane Hydrogen % if cytosolic_ATP f(23) = (-J_ATP - J_ATPK)/W_c; % ATP levels established by mitochondria & glycolysis f(24) = (-J_ADP + J_ATPK)/W_c; % reviewer comment else f(23) = 0; % No ATP change for isolated mitochondria f(24) = 0; % reviewer comment end % %************************************************************************* % Scale Time to minutes (parameters in seconds, but output in minutes) %************************************************************************* % f = 60*f'; % Scale to minutes % %************************************************************************* % Output variables %************************************************************************* % C(1) = J_DH; % Dehydrogenase flux (input function) C(2) = J_C1; % Respiration flux Complex I C(3) = J_C3; % Respiration flux Complex III C(4) = J_C4; % Respiration flux Complex IV C(5) = J_F1; % ATP synthase flux C(6) = J_ANT; % ATP transfer matrix IMS C(7) = J_ATP; % ATP transfer IMS cytosol C(8) = J_ADP; % ADP transfer IMS cytosol C(9) = H_x; % matrix Hydrogen C(10) = Ctot; % Total available cyt-c C(11) = Qtot; % Total available Complex1 C(12) = dPsi-2.3026*(log10(H_x) - pH_e); % unused C(13) = dPsi; % Mitochondrial Membrane potential C(14) = dG_H/F; % Membrame Protomotive Force C(15) = J_Hle; % Proton leak flux C(16) = ATP_e/(ADP_e+ATP_e); % Cytosolic/Medium ATP C(17) = -log10(H_x) - pH_e; % delta matrix and cytosolic pH toth = (4+1)*J_C1 +(4-2)*J_C3 +(2+2)*J_C4 + J_KH ; % total proton extrusion C(18) = 100*((4+1)*J_C1 + (4-2)*J_C3 + (2+2)*J_C4)/toth; % proton extrusion by respiration C(19) = 100*(J_KH)/toth; % proton generation by potassium hydrogen C(20) = 100*(J_DH)/toth; % proton consumption by inpt flux C(21) = 100*(n_A-1)*J_F1/toth; % proton consumption by respiration C(22) = 100*2*J_Pi1/toth; % proton consumption by phosphate-hydrogen C(23) = 100*J_Hle/toth; % proton consumption by proton leaks C(24) = ATP_e; % cytosolic ATP C(25) = ADP_e; % cytosolic ADP end