is_Kset_locally_stable.m

function stable_or_nah = is_Kset_locally_stable(Network_Data, K, fractions)

%%%%%%%%%%%%%%%%%%%%%% Define Model Characteristics %%%%%%%%%%%%%%%%%%%%%%%

Sr = Network_Data.conserved_model_info.Sr;
Lo = Network_Data.conserved_model_info.Lo;
metab_index_old_to_new = ...
    Network_Data.conserved_model_info.metab_index_new_to_old;
no_independent_metabs = size(Sr,1);
L = Network_Data.conserved_model_info.L;
no_metabs = size(Network_Data.S,1);       

initial_conc = [ones(no_metabs,1); fractions];
new_initial_conc = initial_conc;
new_initial_conc(metab_index_old_to_new) = initial_conc;
    
T = new_initial_conc((no_independent_metabs+1):end) - ...
    Lo*new_initial_conc(1:no_independent_metabs);

elasticity_coeff = str2func(strcat('elasticity_coeff_', ...
    Network_Data.model_size));

%%%%%%%%%%%%%% Calculate Jacobian with sampled parameter set %%%%%%%%%%%%%%
%%%%%%%%%%%%%%% @ reference steady state flux distribution %%%%%%%%%%%%%%%%

df_dx = full(jacobian(0,new_initial_conc(1:no_independent_metabs), ...
    K,Sr,L,T,Lo, elasticity_coeff));
eigen_values = eig(df_dx);                                                  % Calculate eigenvalues of Jacobian
real_parts = real(eigen_values);                                            % Identify real part values of eigenvalues
    
if max(real_parts) <= -0.001                                                % If all real parts are negative, the parameter set is marked as locally stable
    stable_or_nah = 1;
else
    stable_or_nah = 0;
end

end