From e2ee569ab8510b8cc713cbd6fd200813f80d6a4f Mon Sep 17 00:00:00 2001
From: agentmess <peder.larson@ucsf.edu>
Date: Wed, 7 Feb 2018 16:56:37 -0800
Subject: [PATCH 1/6] Removing old fitting functions

---
 kinetic_modeling/fit_kPL.m                 |   4 +-
 kinetic_modeling/fit_kinetics_vfa_sameT1.m | 101 ---------------------
 kinetic_modeling/test_fit_kinetics_fcn.m   |  55 -----------
 3 files changed, 3 insertions(+), 157 deletions(-)
 delete mode 100644 kinetic_modeling/fit_kinetics_vfa_sameT1.m
 delete mode 100644 kinetic_modeling/test_fit_kinetics_fcn.m

diff --git a/kinetic_modeling/fit_kPL.m b/kinetic_modeling/fit_kPL.m
index 31dd08f..cdfb271 100644
--- a/kinetic_modeling/fit_kPL.m
+++ b/kinetic_modeling/fit_kPL.m
@@ -186,7 +186,9 @@
     Sfit = reshape(Sfit, [Nx, Nt]);
     ufit = reshape(ufit, [Nx, Nt]);
     objective_val = reshape(objective_val, Nx);
-    disp('100 % complete')
+	if plot_flag
+    	disp('100 % complete')
+    end
 end
 
 end
diff --git a/kinetic_modeling/fit_kinetics_vfa_sameT1.m b/kinetic_modeling/fit_kinetics_vfa_sameT1.m
deleted file mode 100644
index 656a5c6..0000000
--- a/kinetic_modeling/fit_kinetics_vfa_sameT1.m
+++ /dev/null
@@ -1,101 +0,0 @@
-function [R1fit KPXfit  Sfit Minitfit] = fit_kinetics_vfa_sameT1(S, TR, flips, R1est, KPXest)
-% fit_kinetics_vfa_sameT1 - Kinetic Model fitting routine that accounts for 
-% arbitrary RF flip angles over time and between metabolites.
-% This uses the following assumptions to provide more robust fit:
-%   - No fitting of input function, fitting signal starts after input is
-%   complete
-%   - uni-directional conversion from substrate to 1 or 2 metabolic
-%   products (i.e. pyruvate to lactate and pyruvate to alanine)
-%   - Same T1 relaxation rate for all metabolites
-%
-% [R1fit KPXfit  Sfit Minitfit] = fit_kinetics_vfa_sameT1(S, TR, flips, R1est, KPXest)
-%
-% INPUTS
-%	S - signal dynamics [# of metabolites, # of time points]
-%   TR - repetition time per time point
-%   flips - all flip angles [# of metabolites, # of time points x # of phase encodes]
-%	R1est (optional) - relaxation rate initial guess (1/s)
-%	kPXest (optional) - pyruvate to metabolites conversion rate initial guess (1/s)
-% OUTPUTS
-%   R1fit - fit relaxation rates (1/s)
-%   KPXfit - fit conversion rates (1/s)
-%   Sfit - fit signal dynamics
-%   Minitfit - fit initial magnetizations
-%
-% Author: Peder E. Z. Larson
-%
-% (c)2013-2016 The Regents of the University of California.
-% All Rights Reserved.
-
-Nt = size(S,2);
-Nmets = size(S,1);
-if Nmets < 2 || Nmets > 3
-    error('Only fitting of substrate to 1 or 2 metabolic products supported');
-end
-
-[Sscale Mzscale] = flips_scaling_factors(flips, Nt);
-
-Mz_m0 = S(:,1) ./ Sscale(:,1);
-
-if nargin < 5
-    KPXest = 0.02*ones(1,Nmets-1);
-end
-
-if nargin < 4
-    R1est = 1/25; % 1/s
-end
-
-% fit vector - initial guess
-X0 = [R1est KPXest Mz_m0.']; 
-
-    function Sest = model_exchange(x)
-        
-        if Nmets == 3
-            A = [-x(1)-x(2)-x(3) 0 0;
-                +x(2) -x(1) 0;
-                +x(3) 0 -x(1)];
-        else
-            A = [-x(1)-x(2) 0;
-                +x(2) -x(1)]; 
-        end
-        
-        Mz_m(1:Nmets,1) = x(end-Nmets+1:end) ;
-%        Mz_m(1:2,1) = Mz_m0;
-        for It = 1:Nt*Nflips
-             Mz_p(:,It) = Mz_m(:,It) .* cos(flips(:,It));
-             Mxy(:,It) = Mz_m(:,It) .* sin(flips(:,It));
-             Mz_m(:,It+1) = expm(A * TR/Nflips) * Mz_p(:,It);
-        end
-        Sest = squeeze(sum(reshape(Mxy,Nmets,Nflips,Nt),2));
-%             Mz_p(:,It) = Mz_m(:,It) .* Mzscale(:,It);
-%             Sest(:,It) = Mz_m(:,It) .* Sscale(:,It);
-%             Mz_m(:,It+1) = expm(A * TR) * Mz_p(:,It);
-%             
-%         for It = 1:Nt
-%         end
-    end
-
-    function res = g(x)
-        res = model_exchange(x) - S;
-        res = res(:);
-    end
-
-% some refinement possible here
-opts = optimset('Display', 'off');%, 'FinDiffRelStep', [1 1 1 .1 .1]*sqrt(eps)); 
-% FinDiffRelStep - use for normalization? or to preferentially fit certain parameters
-% ie scale step size especially  for MzM0.  Or else scale data
-% FinDiffType
-%'MaxFunEvals', 1e8, 'TolFun', min(abs(S(:,end)))/1e14,'TolX', 1e-14);
-
-%X = fminsearch(@g, X0, opts);
-% set some bounds on the fitting
-% 1/65 < R1 < 1/10, 0 < kPX < 1, 0 < Minit
-lb = [1/65, zeros(1, Nmets-1), zeros(1,Nmets)];
-ub = [1/10, ones(1, Nmets-1), Inf*ones(1,Nmets)];
-
-Xfit = lsqnonlin(@g, X0, lb, ub, opts);
-R1fit = Xfit(1); KPXfit = Xfit([2:Nmets]);
-Minitfit = Xfit(end-Nmets+1:end);
-Sfit = model_exchange(Xfit);
-
-end
\ No newline at end of file
diff --git a/kinetic_modeling/test_fit_kinetics_fcn.m b/kinetic_modeling/test_fit_kinetics_fcn.m
deleted file mode 100644
index 2dfcb45..0000000
--- a/kinetic_modeling/test_fit_kinetics_fcn.m
+++ /dev/null
@@ -1,55 +0,0 @@
-% Script for testing fit_kinetics_* kinetic model fitting functions
-
-clear all
-
-% Test values
-Tin = 6; Tacq = 30; TR = 3; N = Tacq/TR;
-R1 = [1/25 1/25]; KPL = 0.05; std_noise = 0.01;
-k12 = 0.05; Minit = [1; 0.07];
-
-% Test over multiple combinations of flip angle schemes
-flips(1:2,1:N,1) = ones(2,N)*30*pi/180;  % constant, single-band
-flips(1:2,1:N,2) = repmat([20;35]*pi/180,[1 N]);  % constant, multi-band
-flips(1:2,1:N,3) = repmat(vfa_const_amp(N, pi/2), [2,1]);  % RF compensated variable flip angle
-flips(1:2,1:N,4) = [vfa_const_amp(N, pi/2, exp(-TR * ( k12))); ... % T1-effective variable flip angle
-		    vfa_const_amp(N, pi/2, exp(-TR * ( - k12)))];
-flips(1:2,1:N,5) = [vfa_const_amp(N, pi/2, exp(-TR * ( k12))); ... % max lactate SNR variable flip angle
-		    vfa_opt_signal(N, exp(-TR * ( R1(2))))];
-flips(1:2,1:N,6) = [vfa_const_amp(N, pi/2, exp(-TR * (k12))); ... % saturation recovery
-		    ones(1,N)*pi/2];
-
-         t = [0:N-1]*TR + Tin;
-
-        figure
- subplot(121) , plot(t, squeeze(flips(1,:,:))*180/pi)
-title('Pyruvate flips')
- subplot(122) , plot(t, squeeze(flips(2,:,:))*180/pi)
-title('Lactate flips')
-legend('constant','multiband', 'vfa', 'T1-effective vfa', 'max lactate SNR vfa', 'Saturation Recovery')
-        
-%% Test fitting
- for Iflips = 1:size(flips,3)
-     [Mxy Mz] = simulate_2site_model(Tin, R1, [KPL 0], flips(:,:,Iflips), TR);
-     % no noise
-    [R1fit(Iflips) KPLfit(Iflips) Sfit(1:size(Mxy,1), 1:size(Mxy,2),  Iflips)] = fit_kinetics_vfa_sameT1(Mxy, TR, flips(:,:,Iflips));
-    % add noise
-    Sn(1:size(Mxy,1), 1:size(Mxy,2),  Iflips) = Mxy + randn(size(Mxy))*std_noise;
-    [R1fitn(Iflips) KPLfitn(Iflips) Snfit(1:size(Mxy,1), 1:size(Mxy,2),  Iflips)] = fit_kinetics_vfa_sameT1(Sn(:,:,Iflips), TR, flips(:,:,Iflips));
- end
-
-disp(sprintf('Input R1 = %f (pyr) %f (lac), kPL = %f', R1(1), R1(2), KPL))
-disp('Noiseless fit results:')
-disp(['R1fit  = ' num2str(R1fit)])
-disp(['KPLfit = ' num2str(KPLfit)])
-disp('Noisy fit results:')
-disp(['R1fitn  = ' num2str(R1fitn)])
-disp(['KPLfitn = ' num2str(KPLfitn)])
- 
-  figure
- subplot(121) , plot(t, squeeze(Sn(1,:,:)))
-    hold on, plot(t, squeeze(Snfit(1,:,:)),'--')
-title('Pyruvate signals and fits (dashed)')
- subplot(122) , plot(t, squeeze(Sn(2,:,:)))
-    hold on, plot(t, squeeze(Snfit(2,:,:)),'--')
-title('Lactate signals and fits (dashed)')
-legend('constant','multiband', 'vfa', 'T1-effective vfa', 'max lactate SNR vfa', 'Saturation Recovery')

From 94f5afa5cb63907b03b0f4ea19b6d79017188e7c Mon Sep 17 00:00:00 2001
From: agentmess <peder.larson@ucsf.edu>
Date: Tue, 13 Feb 2018 10:27:10 -0800
Subject: [PATCH 2/6] Initial commit of new functions

Not complete or tested
---
 kinetic_modeling/fit_HP_kinetics.m | 226 +++++++++++++++++++++++++++++
 kinetic_modeling/trajectories.m    |  49 +++++++
 2 files changed, 275 insertions(+)
 create mode 100644 kinetic_modeling/fit_HP_kinetics.m
 create mode 100644 kinetic_modeling/trajectories.m

diff --git a/kinetic_modeling/fit_HP_kinetics.m b/kinetic_modeling/fit_HP_kinetics.m
new file mode 100644
index 0000000..394719d
--- /dev/null
+++ b/kinetic_modeling/fit_HP_kinetics.m
@@ -0,0 +1,226 @@
+function [params_fit, Sfit, ufit, objective_val] = fit_HP_kinetics(S, TR, flips, model, params_fixed, params_est, noise_level, plot_flag)
+% fit_HP_kinetics - Kinetic model fitting function for HP 13C MRI.
+%
+% Models - 'inputless' (default), 'boxcar-input', 'gamma-input'
+% 13C of conversion rate by fitting of
+% product (e.g. lactate) signal at each time point.  Substrate (e.g.
+% pyruvate) signal taken as is, and not fit to any function, eliminating
+% need to make any assumptions about the input function.
+% This uses the following assumptions:
+%   - uni-directional conversion from substrate to metabolic products (i.e.
+%   pyruvate to lactate)
+%   - initial lactate magnetization is zero (need to add)
+% It also allows for fixing of parameters. Based on simulations, our
+% current recommendation is to fix pyruvate T1, as it doesn't impact kPL substantially.
+%
+% [params_fit, Sfit, ufit, objective_val] = fit_kPL(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
+%
+% All params_* values are structures, with possible fields of 'kPL', 'R1L',
+% and 'R1P', and units of 1/s.
+% INPUTS
+%	S - signal dynamics [voxels, # of metabolites, # of time points]
+%   TR - repetition time per time point flips - all flip angles [# of
+%   metabolites, # of time points x # of phase encodes]
+%	params_fixed - structure of fixed parameters and values (1/s).  parameters not in
+%       this structure will be fit
+%   params_est (optional) - structure of estimated values for fit parameters pyruvate to metabolites conversion rate initial guess (1/s)
+%       Also can include upper and lower bounds on parameters as *_lb and
+%       *_ub (e.g. R1L_lb, R1L_ub)
+%   noise_level (optional) - estimate standard deviation of noise in data
+%       to use maximum likelihood fit of magnitude data (with Rician noise
+%       distribution)
+%   plot_flag (optional) - plot fits
+% OUTPUTS
+%   params_fit - structure of fit parameters 
+%   Sfit - fit curves
+%   ufit - derived input function (unitless)
+%   objective_val - measure of fit error
+%
+% EXAMPLES - see test_fit_HP_kinetics.m
+%
+% Authors: John Maidens,  Peder E. Z. Larson
+%
+% (c)2015-2018 The Regents of the University of California. All Rights
+% Reserved.
+
+if nargin < 4 || isempty(model)
+    model = 'inputless';
+end
+
+params_all = {'kPL', 'R1L', 'R1P', 'L0_start', 'Rinj', 'Tarrival', 'Tbolus'};
+params_default_est = [0.02, 1/25, 1/25, 0 0.1 0 8];
+params_default_lb = [0, 1/60, 1/60, 0 0 -30 0];
+params_default_ub = [Inf, 1/10, 1/10, Inf Inf 30 Inf];
+
+if nargin < 5 || isempty(params_fixed)
+    params_fixed = struct([]);
+end
+
+if nargin < 6 || isempty(params_est)
+    params_est = struct([]);
+end
+
+I_params_est = [];
+for n = 1:length(params_all)
+    if ~isfield(params_fixed, params_all(n))
+        I_params_est = [I_params_est, n];
+    end
+end
+Nparams_to_fit = length(I_params_est);
+
+for n = 1:Nparams_to_fit
+    param_name = params_all{I_params_est(n)};
+    if isfield(params_est, param_name)
+        params_est_vec(n) = params_est.(param_name);
+    else
+        params_est_vec(n) = params_default_est(I_params_est(n));
+    end
+    if isfield(params_est, [param_name '_lb'])
+        params_lb(n) = params_est.([param_name '_lb']);
+    else
+        params_lb(n) = params_default_lb(I_params_est(n));
+    end
+    if isfield(params_est, [param_name '_ub'])
+        params_ub(n) = params_est.([param_name '_ub']);
+    else
+        params_ub(n) = params_default_ub(I_params_est(n));
+    end
+end
+
+
+if nargin < 6 || isempty(noise_level)
+    % no noise level provided, so use least-squares fit (best for Gaussian
+    % zero-mean noise)
+    fit_method = 'ls';
+else
+    % otherwise use maximum likelihood (good for Rician noise from
+    % magnitudes)
+    fit_method = 'ml';
+end
+
+if nargin < 7
+    plot_flag = 0;
+end
+
+if plot_flag
+    disp('==== Computing parameter map ====')
+end
+
+size_S = size(S);  ndimsx = length(size_S)-2;
+Nt = size_S(end); t = [0:Nt-1]*TR;
+Nx = size_S(1:ndimsx);
+if isempty(Nx)
+    Nx = 1;
+end
+S = reshape(S, [prod(Nx), 2, Nt]);  % put all spatial locations in first dimension
+
+[Sscale, Mzscale] = flips_scaling_factors(flips, Nt);
+
+params_fit_vec = zeros([prod(Nx),Nparams_to_fit]);  objective_val = zeros([1,prod(Nx)]);
+Sfit = zeros([prod(Nx),Nt]); ufit = zeros([prod(Nx),Nt]);
+
+for i=1:size(S, 1)
+    if length(Nx) > 1 && plot_flag
+        disp([num2str( floor(100*(i-1)/size(S, 1)) ) '% complete'])
+    end
+    % observed magnetization (Mxy)
+    y1 = reshape(S(i, 1, :), [1, Nt]); % pyr
+    y2 = reshape(S(i, 2, :), [1, Nt]); % lac
+    if any(y1 ~= 0)
+        % % plot of observed data for debugging
+        % figure(1)
+        % plot(t, y1, t, y2)
+        % xlabel('time (s)')
+        % ylabel('measured magnetization (au)')
+        % legend('pyruvate', 'lactate')
+        
+        % estimate state magnetization (MZ) based on scaling from RF pulses
+        x1 = y1./Sscale(1, :);
+        x2 = y2./Sscale(2, :);
+        
+        % fit to data
+        options = optimoptions(@fminunc,'Display','none','Algorithm','quasi-newton');
+        lsq_opts = optimset('Display','none','MaxIter', 500, 'MaxFunEvals', 500);
+        switch(fit_method)
+            case 'ls'
+                obj = @(var) (x2 - trajectories_frompyr(var, x1, Mzscale, params_fixed, TR)) .* Sscale(2,:);  % perform least-squares in signal domain
+                [params_fit_vec(i,:),objective_val(i)] = lsqnonlin(obj, params_est_vec, params_lb, params_ub, lsq_opts);
+                
+            case 'ml'
+                obj = @(var) negative_log_likelihood_rician_frompyr(var, x1, x2, Mzscale, params_fixed, TR, noise_level.*(Sscale(2,:).^2));
+                [params_fit_vec(i,:), objective_val(i)] = fminunc(obj, params_est_vec, options);
+                
+        end
+        [Sfit(i,:), ufit(i,:)] = trajectories_frompyr(params_fit_vec(i,:), x1, Mzscale, params_fixed, TR);
+        Sfit(i,:) = Sfit(i,:)  .* Sscale(2, :);
+        ufit(i,:) = ufit(i,:)  .* Sscale(1, :);
+        
+        if plot_flag
+            % plot of fit for debugging
+            figure(99)
+            subplot(2,1,1)
+            plot(t, x1, t, x2, t, Sfit(i,:)./ Sscale(2, :),'--', t, ufit(i,:)./ Sscale(1, :), 'k:')
+            xlabel('time (s)')
+            ylabel('state magnetization (au)')
+            subplot(2,1,2)
+            plot(t, y1, t, y2, t, Sfit(i,:),'--', t, ufit(i,:), 'k:')
+            xlabel('time (s)')
+            ylabel('signal (au)')
+            title(num2str(params_fit_vec(i,1:end-1),4)) % don't display L0_start value
+            legend('pyruvate', 'lactate', 'lactate fit', 'input estimate')
+            drawnow, pause(0.5)
+        end
+    end
+end
+
+
+params_fit = struct([]);
+nfit = 0;
+for n = 1:length(params_all)-1  % don't output L0_start
+    if ~isfield(params_fixed, params_all(n))
+        nfit = nfit+1;
+        params_fit(1).(params_all{n})= params_fit_vec(:,nfit);
+    end
+end
+
+if length(Nx) > 1
+    for n = 1:Nparams_to_fit-1 % don't output L0_start
+        param_name = params_all{I_params_est(n)};
+        params_fit.(param_name) = reshape(params_fit.(param_name), Nx);
+    end
+    
+    
+    Sfit = reshape(Sfit, [Nx, Nt]);
+    ufit = reshape(ufit, [Nx, Nt]);
+    objective_val = reshape(objective_val, Nx);
+	if plot_flag
+    	disp('100 % complete')
+    end
+end
+
+end
+
+function [ l1 ] = negative_log_likelihood_rician_frompyr(params_fit, x1, x2, Mzscale, params_fixed, TR, noise_level)
+%FUNCTION NEGATIVE_LOG_LIKELIHOOD_RICIAN Computes log likelihood for
+%    compartmental model with Rician noise
+% noise level is scaled for state magnetization (Mz) domain
+
+N = size(x1,2);
+
+% compute trajectory of the model with parameter values
+x2fit = trajectories_frompyr(params_fit, x1, Mzscale, params_fixed, TR);
+
+% compute negative log likelihood
+l1 = 0;
+for t = 1:N
+    for k = 1
+        l1 = l1 - (...
+            log(x2(k, t)) - log(noise_level(t)) ...
+            - (x2(k, t)^2 + x2fit(k, t)^2)/(2*noise_level(t)) ...
+            + x2(k, t)*x2fit(k, t)/noise_level(t) ...
+            + log(besseli(0, x2(k, t)*x2fit(k, t)/noise_level(t), 1))...
+            );
+    end
+end
+end
+
diff --git a/kinetic_modeling/trajectories.m b/kinetic_modeling/trajectories.m
new file mode 100644
index 0000000..77a9fb1
--- /dev/null
+++ b/kinetic_modeling/trajectories.m
@@ -0,0 +1,49 @@
+function [x1, x2] = trajectories_withgammainput( params_fit, params_fixed , TR, N, Mzscale )
+
+% all using MZ component so can account for variable flip angles
+
+x1 = zeros(1, N); x2 = zeros(1,N);
+
+params_all = {'kPL', 'R1L', 'R1P', 'Rinj', 'Tarrival', 'A','B'};
+nfit = 0;
+for n = 1:length(params_all)
+    if isfield(params_fixed, params_all(n))
+        eval([params_all{n} '= params_fixed.(params_all{n});']);
+    else
+        nfit = nfit+1;
+        eval([params_all{n} '= params_fit(nfit);']);
+    end
+end
+
+for It=1:N
+    
+    
+    t = (It-1)*TR;  % solving for Mz at time t (accounting for previous TR interval)
+    
+    if It == 1
+        P0 = 0;
+        L0 = 0;
+        if Tarrival < 0
+            % account for longer period of bolus prior to imaging
+            t_preimaging = [floor(Tarrival/TR):0]*TR; % t = 0
+           u = sum( gampdf(t_preimaging-Tarrival,A,B)*Rinj );
+        else
+           u = gampdf(t-Tarrival,A,B)*Rinj;
+        end
+        
+    else
+        P0 = x1(It-1)*Mzscale(1, It-1);
+        L0 = x2(It-1)*Mzscale(2, It-1);
+        u = gampdf(t-Tarrival,A,B)*Rinj;
+    end
+        
+    % solve next time point under assumption of constant input during TR
+    x2(It) = exp(-R1L*TR)*((L0*R1L*R1P - L0*R1L^2 - kPL*u + L0*R1L*kPL + P0*R1L*kPL)/(R1L*(R1P - R1L + kPL)) + (kPL*u*exp(R1L*TR))/(R1L*(R1P - R1L + kPL))) - exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)));
+    
+    % solution for next source (pyruvate) time point if needed:
+    x1(It) = (exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)))*(R1P - R1L + kPL))/kPL;
+    
+end
+
+end
+

From d8e90b0ca54bebda2d3f5385567e2a08455192a7 Mon Sep 17 00:00:00 2001
From: agentmess <peder.larson@ucsf.edu>
Date: Tue, 20 Mar 2018 08:55:38 -0700
Subject: [PATCH 3/6] checkpoint commit, still lots to do

---
 .../{fit_HP_kinetics.m => fit_pyr_kinetics.m} | 118 +++++++++++++-----
 kinetic_modeling/trajectories.m               |  62 ++++++---
 .../demo_dynamic_mrs_multishot_fit.m          |  19 ++-
 3 files changed, 151 insertions(+), 48 deletions(-)
 rename kinetic_modeling/{fit_HP_kinetics.m => fit_pyr_kinetics.m} (66%)

diff --git a/kinetic_modeling/fit_HP_kinetics.m b/kinetic_modeling/fit_pyr_kinetics.m
similarity index 66%
rename from kinetic_modeling/fit_HP_kinetics.m
rename to kinetic_modeling/fit_pyr_kinetics.m
index 394719d..5f08063 100644
--- a/kinetic_modeling/fit_HP_kinetics.m
+++ b/kinetic_modeling/fit_pyr_kinetics.m
@@ -1,4 +1,4 @@
-function [params_fit, Sfit, ufit, objective_val] = fit_HP_kinetics(S, TR, flips, model, params_fixed, params_est, noise_level, plot_flag)
+function [params_fit, Sfit, ufit, objective_val] = fit_pyr_kinetics(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
 % fit_HP_kinetics - Kinetic model fitting function for HP 13C MRI.
 %
 % Models - 'inputless' (default), 'boxcar-input', 'gamma-input'
@@ -15,10 +15,11 @@
 %
 % [params_fit, Sfit, ufit, objective_val] = fit_kPL(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
 %
-% All params_* values are structures, with possible fields of 'kPL', 'R1L',
-% and 'R1P', and units of 1/s.
+% All params_* values are structures, with possible fields of 'kPL' (1/s), 'R1L' (1/s),
+% and 'R1P', and units of 1/s.  MORE
 % INPUTS
 %	S - signal dynamics [voxels, # of metabolites, # of time points]
+%		Substrate (e.g. Pyruvate) should be the first metabolite, followed by each product
 %   TR - repetition time per time point flips - all flip angles [# of
 %   metabolites, # of time points x # of phase encodes]
 %	params_fixed - structure of fixed parameters and values (1/s).  parameters not in
@@ -47,10 +48,30 @@
     model = 'inputless';
 end
 
-params_all = {'kPL', 'R1L', 'R1P', 'L0_start', 'Rinj', 'Tarrival', 'Tbolus'};
-params_default_est = [0.02, 1/25, 1/25, 0 0.1 0 8];
-params_default_lb = [0, 1/60, 1/60, 0 0 -30 0];
-params_default_ub = [Inf, 1/10, 1/10, Inf Inf 30 Inf];
+size_S = size(S);  ndimsx = length(size_S)-2;
+Nt = size_S(end); t = [0:Nt-1]*TR;
+Nx = size_S(1:ndimsx);
+Nmets = size_S(end-1);
+if isempty(Nx)
+    Nx = 1;
+end
+
+params_all = {'kPL', 'kPB', 'kPA', ...
+	'R1P', 'R1L', 'R1A', 'R1B', ...
+	'S0_L', 'S0_B', 'S0_A', ...
+	'Rinj', 'Tarrival', 'Tbolus'};
+params_default_est = [0.01, 0.01, 0.01, ...
+	1/30, 1/25, 1/25, 1/15, ...
+	0, 0, 0, ...
+	0.1, 0, 8];
+params_default_lb = [-Inf, -Inf, -Inf, ...
+	 1/50, 1/50, 1/50, 1/50, ...
+	 -Inf, -Inf, -Inf, ...
+	  0, -30, 0];
+params_default_ub = [Inf, Inf, Inf, ...
+	1/10, 1/10, 1/10, 1/5 , ...
+	 Inf, Inf, Inf, ...
+	 Inf 30 Inf];
 
 if nargin < 5 || isempty(params_fixed)
     params_fixed = struct([]);
@@ -60,6 +81,16 @@
     params_est = struct([]);
 end
 
+% Supports up to 3 metabolic products (e.g. alanine, lactate, bicarb)
+switch Nmets
+	case 2 % assume pyruvate & lactate
+params_fixed.kPA = 0;  params_fixed.S0_A = 0;
+params_fixed.kPB = 0;  params_fixed.S0_B = 0;
+	case 3 % assume pyruvate & lactate & bicarbonate
+params_fixed.kPA = 0;   params_fixed.S0_A = 0;
+end
+
+
 I_params_est = [];
 for n = 1:length(params_all)
     if ~isfield(params_fixed, params_all(n))
@@ -106,33 +137,21 @@
     disp('==== Computing parameter map ====')
 end
 
-size_S = size(S);  ndimsx = length(size_S)-2;
-Nt = size_S(end); t = [0:Nt-1]*TR;
-Nx = size_S(1:ndimsx);
-if isempty(Nx)
-    Nx = 1;
-end
-S = reshape(S, [prod(Nx), 2, Nt]);  % put all spatial locations in first dimension
+S = reshape(S, [prod(Nx), Nmets, Nt]);  % put all spatial locations in first dimension
 
 [Sscale, Mzscale] = flips_scaling_factors(flips, Nt);
 
 params_fit_vec = zeros([prod(Nx),Nparams_to_fit]);  objective_val = zeros([1,prod(Nx)]);
-Sfit = zeros([prod(Nx),Nt]); ufit = zeros([prod(Nx),Nt]);
+Sfit = zeros([prod(Nx),Nmets,Nt]); ufit = zeros([prod(Nx),Nt]);
 
 for i=1:size(S, 1)
     if length(Nx) > 1 && plot_flag
         disp([num2str( floor(100*(i-1)/size(S, 1)) ) '% complete'])
     end
     % observed magnetization (Mxy)
-    y1 = reshape(S(i, 1, :), [1, Nt]); % pyr
-    y2 = reshape(S(i, 2, :), [1, Nt]); % lac
+    Mxy = reshape(S(i, :, :), [Nmets, Nt]); % pyr
+
     if any(y1 ~= 0)
-        % % plot of observed data for debugging
-        % figure(1)
-        % plot(t, y1, t, y2)
-        % xlabel('time (s)')
-        % ylabel('measured magnetization (au)')
-        % legend('pyruvate', 'lactate')
         
         % estimate state magnetization (MZ) based on scaling from RF pulses
         x1 = y1./Sscale(1, :);
@@ -143,15 +162,16 @@
         lsq_opts = optimset('Display','none','MaxIter', 500, 'MaxFunEvals', 500);
         switch(fit_method)
             case 'ls'
-                obj = @(var) (x2 - trajectories_frompyr(var, x1, Mzscale, params_fixed, TR)) .* Sscale(2,:);  % perform least-squares in signal domain
+                obj = @(var) (x2 - trajectories_inputless(model, var, x1, Mzscale, params_fixed, TR)) .* Sscale(2,:);  % perform least-squares in signal domain
                 [params_fit_vec(i,:),objective_val(i)] = lsqnonlin(obj, params_est_vec, params_lb, params_ub, lsq_opts);
                 
             case 'ml'
-                obj = @(var) negative_log_likelihood_rician_frompyr(var, x1, x2, Mzscale, params_fixed, TR, noise_level.*(Sscale(2,:).^2));
+                obj = @(var) negative_log_likelihood_rician_inputless(model, var, x1, x2, Mzscale, params_fixed, TR, noise_level.*(Sscale(2,:).^2));
                 [params_fit_vec(i,:), objective_val(i)] = fminunc(obj, params_est_vec, options);
                 
         end
-        [Sfit(i,:), ufit(i,:)] = trajectories_frompyr(params_fit_vec(i,:), x1, Mzscale, params_fixed, TR);
+        
+        [Sfit(i,:), ufit(i,:)] = trajectories(params_fit_vec(i,:), x1, Mzscale, params_fixed, TR);
         Sfit(i,:) = Sfit(i,:)  .* Sscale(2, :);
         ufit(i,:) = ufit(i,:)  .* Sscale(1, :);
         
@@ -200,7 +220,7 @@
 
 end
 
-function [ l1 ] = negative_log_likelihood_rician_frompyr(params_fit, x1, x2, Mzscale, params_fixed, TR, noise_level)
+function [ l1 ] = negative_log_likelihood_rician_inputless(params_fit, x1, x2, Mzscale, params_fixed, TR, noise_level)
 %FUNCTION NEGATIVE_LOG_LIKELIHOOD_RICIAN Computes log likelihood for
 %    compartmental model with Rician noise
 % noise level is scaled for state magnetization (Mz) domain
@@ -208,7 +228,7 @@
 N = size(x1,2);
 
 % compute trajectory of the model with parameter values
-x2fit = trajectories_frompyr(params_fit, x1, Mzscale, params_fixed, TR);
+x2fit = trajectories_inputless(params_fit, x1, Mzscale, params_fixed, TR);
 
 % compute negative log likelihood
 l1 = 0;
@@ -224,3 +244,45 @@
 end
 end
 
+function [Mz_products, u] = trajectories_inputless( params_fit, Mz_pyr, Mzscale, params_fixed , TR )
+% Compute product magnetizations using a uni-directional two-site model
+% Uses substrate magnetization measurements, estimated relaxation and
+% conversion rates
+
+Nmets = size(Mzscale,1); N = size(Mzscale,2);
+Mz_all = zeros(Nmets, N);
+u = zeros(1,N);
+
+params_all = {'kPL', 'kPB', 'kPA', ...
+	'R1P', 'R1L', 'R1A', 'R1B', ...
+	'S0_L', 'S0_B', 'S0_A'};
+
+nfit = 0;
+for n = 1:length(params_all)
+    if isfield(params_fixed, params_all(n))
+        eval([params_all{n} '= params_fixed.(params_all{n});']);
+    else
+        nfit = nfit+1;
+        eval([params_all{n} '= params_fit(nfit);']);
+    end
+end
+
+Mz_all(1,:) = Mz_pyr;
+Mz_all(2,1) = S0_L;
+Mz_all(3,1) = S0_B;
+Mz_all(4,1) = S0_A;
+
+for It=1:N-1
+    
+    Mz_init = Mz_all(:,It) .* Mzscale(:, It);
+    
+    % estimate input, assuming this is constant during TR interval
+    u(It) = ( Mz_pyr(It+1) - Mz_init(1)*exp((- R1P - kPL - kPB - kPA)*TR) ) * (R1P + kPL) / (1 - exp((- R1P - kPL - kPB - kPA)*TR));
+    
+    % solve next time point under assumption of constant input during TR
+    x2(It+1) = exp(-R1L*TR)*((L0*R1L*R1P - L0*R1L^2 - kPL*u(It) + L0*R1L*kPL + P0*R1L*kPL)/(R1L*(R1P - R1L + kPL)) + (kPL*u(It)*exp(R1L*TR))/(R1L*(R1P - R1L + kPL))) - exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u(It) + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u(It)*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)));
+    
+    
+end
+
+end
diff --git a/kinetic_modeling/trajectories.m b/kinetic_modeling/trajectories.m
index 77a9fb1..37f3dd1 100644
--- a/kinetic_modeling/trajectories.m
+++ b/kinetic_modeling/trajectories.m
@@ -1,10 +1,15 @@
-function [x1, x2] = trajectories_withgammainput( params_fit, params_fixed , TR, N, Mzscale )
+function [Mz u] = trajectories( model, params_fit, params_fixed , TR, Mzscale )
 
 % all using MZ component so can account for variable flip angles
+Nmets = size(Mzscale,1); N = size(Mzscale,2);
+Mz = zeros(Nmets, N);
 
-x1 = zeros(1, N); x2 = zeros(1,N);
 
-params_all = {'kPL', 'R1L', 'R1P', 'Rinj', 'Tarrival', 'A','B'};
+params_all = {'kPL', 'kPB', 'kPA', ...
+	'R1P', 'R1L', 'R1A', 'R1B', ...
+	'S0_L', 'S0_B', 'S0_A', ...
+	'Rinj', 'Tarrival', 'Tbolus'};
+
 nfit = 0;
 for n = 1:length(params_all)
     if isfield(params_fixed, params_all(n))
@@ -15,34 +20,59 @@
     end
 end
 
+% gamma input parameters
+A = Tbolus;
+B = Tbolus/4; % assumption that leads to reasonable looking shape
+
+% evolution matrix
+switch Nmets
+    case 2
+        Amat = [-R1P - kPL, 0; ...
+            kPL, -R1L];
+    case 3
+        Amat = [-R1P - kPL - kPB, 0, 0; ...
+            kPL, -R1L, 0; ...
+            kPB, 0, -R1B];
+    case 4
+        Amat = [-R1P - kPL - kPB - kPA, 0, 0, 0; ...
+            kPL, -R1L, 0, 0; ...
+            kPB, 0, -R1B, 0; ...
+            kPA, 0, 0, -R1A];
+end
+
 for It=1:N
     
-    
+    if strcmp(model, 'inputless')
+        Mz_init = Mz(:,It-1) .* Mzscale(:, It-1);
+        % FIX:
+         u(It) = ( x1(t+1) - P0*exp((- R1P - kPL)*TR) ) * (R1P + kPL) / (1 - exp((- R1P - kPL)*TR));
+    else
+        
     t = (It-1)*TR;  % solving for Mz at time t (accounting for previous TR interval)
     
     if It == 1
-        P0 = 0;
-        L0 = 0;
+    	% initial Mz in TR after RF pulse
+        Mz_init = zeros(NMets,1);
+
         if Tarrival < 0
             % account for longer period of bolus prior to imaging
             t_preimaging = [floor(Tarrival/TR):0]*TR; % t = 0
-           u = sum( gampdf(t_preimaging-Tarrival,A,B)*Rinj );
+           u(It) = sum( gampdf(t_preimaging-Tarrival,A,B)*Rinj );
         else
-           u = gampdf(t-Tarrival,A,B)*Rinj;
+           u(It) = gampdf(t-Tarrival,A,B)*Rinj;
         end
         
     else
-        P0 = x1(It-1)*Mzscale(1, It-1);
-        L0 = x2(It-1)*Mzscale(2, It-1);
-        u = gampdf(t-Tarrival,A,B)*Rinj;
+    	% initial Mz in TR after RF pulse
+    	Mz_init = Mz(:,It-1) .* Mzscale(:, It-1);
+
+        u(It) = gampdf(t-Tarrival,A,B)*Rinj;
+    end
     end
         
     % solve next time point under assumption of constant input during TR
-    x2(It) = exp(-R1L*TR)*((L0*R1L*R1P - L0*R1L^2 - kPL*u + L0*R1L*kPL + P0*R1L*kPL)/(R1L*(R1P - R1L + kPL)) + (kPL*u*exp(R1L*TR))/(R1L*(R1P - R1L + kPL))) - exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)));
-    
-    % solution for next source (pyruvate) time point if needed:
-    x1(It) = (exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)))*(R1P - R1L + kPL))/kPL;
-    
+	Mz(:,It) = expm(Amat*TR)*Mz_init + u; % WRONG WRONG
+	    
 end
 
 end
diff --git a/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m b/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m
index 945cd02..70cb78f 100644
--- a/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m	
+++ b/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m	
@@ -10,6 +10,8 @@
 clear all
 plot_flag = 1;
 int_range = 25;
+    R1P_est = 1/30; R1L_est = 1/25;  R1B_est = 1/15; R1A_est = 1/25;
+    kPL_est = 0.01; kPB_est = .01; kPA_est = .01;
 
 for ratnum = 1:4;
     
@@ -54,18 +56,27 @@
     end
     
     
-    clear params_fixed params_est params_fit kPL
-    R1P_est = 1/30; R1L_est = 1/25;  kPL_est = 0.01;
+    clear params_fixed params_est params_fit kPL kPB kPA
     params_fixed.R1P = R1P_est; params_fixed.R1L = R1L_est;
-    params_est.kPL = kPL_est;
+    params_fixed.R1B = R1B_est; params_fixed.R1A = R1A_est;
+    params_est.kPL = kPL_est; params_est.kPB = kPB_est; params_est.kPA = kPA_est;
     
     for n = 1:Nshots
         load(sprintf('rat%d/shot%d',ratnum,n));
         pyr = real(sum(spectra_dynamic(pyr_center+[-int_range:int_range],:),1));
         lac = real(sum(spectra_dynamic(lac_center+[-int_range:int_range],:),1));
+        ala = real(sum(spectra_dynamic(ala+[-int_range:int_range],:),1));
+        bicarb = real(sum(spectra_dynamic(bicarb+[-int_range:int_range],:),1));
         [params_fit Sfit] = fit_kPL([pyr;lac] , TR, repmat(flip, [2 length(pyr)]), params_fixed, params_est, [], plot_flag);
         kPL(n) = params_fit.kPL;
+        
+        [params_fit Sfit] = fit_pyr_kinetics([pyr;lac;bicarb;ala] , TR, repmat(flip, [4 length(pyr)]), params_fixed, params_est, [], plot_flag);
+        kPL(n) = params_fit.kPL; kPB(n) = params_fit.kPB; kPA(n) = params_fit.kPA;
+        
     end
     
-    disp(['Rat #' int2str(ratnum) ' kPL fits: ' num2str(kPL) ' 1/s'])
+    disp(['Rat #' int2str(ratnum)]);
+    disp(['   kPL fits: ' num2str(kPL) ' 1/s'])
+    disp(['   kPB fits: ' num2str(kPB) ' 1/s'])
+    disp(['   kPA fits: ' num2str(kPA) ' 1/s'])
 end
\ No newline at end of file

From 68dbc9b8fceaffcca40fe3e3a82514ce80217d2f Mon Sep 17 00:00:00 2001
From: agentmess <peder.larson@ucsf.edu>
Date: Tue, 20 Mar 2018 22:23:12 -0700
Subject: [PATCH 4/6] multiple metabolites now supported in testing and with
 data

rat kidneys dynamic data shows good agreement with prior functions
---
 kinetic_modeling/fit_kPL_withgammainput.m     |   2 +-
 kinetic_modeling/fit_pyr_kinetics.m           | 148 +++++++------
 kinetic_modeling/test_fit_pyr_kinetics.m      | 196 ++++++++++++++++++
 .../demo_dynamic_mrs_multishot_fit.m          |  29 +--
 simulations/simulate_Nsite_model.m            |  80 +++++++
 5 files changed, 374 insertions(+), 81 deletions(-)
 create mode 100644 kinetic_modeling/test_fit_pyr_kinetics.m
 create mode 100644 simulations/simulate_Nsite_model.m

diff --git a/kinetic_modeling/fit_kPL_withgammainput.m b/kinetic_modeling/fit_kPL_withgammainput.m
index efc1957..493cbd6 100644
--- a/kinetic_modeling/fit_kPL_withgammainput.m
+++ b/kinetic_modeling/fit_kPL_withgammainput.m
@@ -147,7 +147,7 @@
                 end
                 
                 % Fit all data
-                obj = @(var) trajectory_difference_all(var, x1, x2, params_fixed, TR, Mzscale);
+                obj = @(var) trajectory_difference_all(var, x1, x2, params_fixed, TR, Mzscale);  % add Sscale?
                 [params_fit_vec(i,:),objective_val(i)] = lsqnonlin(obj, params_est_vec2, params_lb, params_ub, lsq_opts);
                 
             case 'ml'
diff --git a/kinetic_modeling/fit_pyr_kinetics.m b/kinetic_modeling/fit_pyr_kinetics.m
index 5f08063..51136ed 100644
--- a/kinetic_modeling/fit_pyr_kinetics.m
+++ b/kinetic_modeling/fit_pyr_kinetics.m
@@ -1,27 +1,25 @@
 function [params_fit, Sfit, ufit, objective_val] = fit_pyr_kinetics(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
-% fit_HP_kinetics - Kinetic model fitting function for HP 13C MRI.
+% fit_pyr_kinetics - Kinetic model fitting function for HP 13C MRI.
 %
-% Models - 'inputless' (default), 'boxcar-input', 'gamma-input'
-% 13C of conversion rate by fitting of
-% product (e.g. lactate) signal at each time point.  Substrate (e.g.
-% pyruvate) signal taken as is, and not fit to any function, eliminating
+% Fits product signals, assuming origination from a single substrate
+% In other words, pyruvate to lactate, bicarbonate and alanine.
+% An input-less method is used, eliminating
 % need to make any assumptions about the input function.
 % This uses the following assumptions:
 %   - uni-directional conversion from substrate to metabolic products (i.e.
 %   pyruvate to lactate)
-%   - initial lactate magnetization is zero (need to add)
 % It also allows for fixing of parameters. Based on simulations, our
-% current recommendation is to fix pyruvate T1, as it doesn't impact kPL substantially.
+% current recommendation is to fix pyruvate T1, as it doesn't impact kPX substantially.
 %
-% [params_fit, Sfit, ufit, objective_val] = fit_kPL(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
+% [params_fit, Sfit, ufit, objective_val] = fit_pyr_kinetics(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
 %
-% All params_* values are structures, with possible fields of 'kPL' (1/s), 'R1L' (1/s),
-% and 'R1P', and units of 1/s.  MORE
+% All params_* values are structures, including possible fields of 'kPL', 'kPB', 'kPA', (1/s),
+% 'R1P', 'R1L', 'R1B', 'R1A' (1/s).
 % INPUTS
 %	S - signal dynamics [voxels, # of metabolites, # of time points]
 %		Substrate (e.g. Pyruvate) should be the first metabolite, followed by each product
-%   TR - repetition time per time point flips - all flip angles [# of
-%   metabolites, # of time points x # of phase encodes]
+%   TR - repetition time per time point
+%	flips - all flip angles [# of metabolites, # of time points x # of phase encodes]
 %	params_fixed - structure of fixed parameters and values (1/s).  parameters not in
 %       this structure will be fit
 %   params_est (optional) - structure of estimated values for fit parameters pyruvate to metabolites conversion rate initial guess (1/s)
@@ -32,7 +30,7 @@
 %       distribution)
 %   plot_flag (optional) - plot fits
 % OUTPUTS
-%   params_fit - structure of fit parameters 
+%   params_fit - structure of fit parameters
 %   Sfit - fit curves
 %   ufit - derived input function (unitless)
 %   objective_val - measure of fit error
@@ -44,9 +42,6 @@
 % (c)2015-2018 The Regents of the University of California. All Rights
 % Reserved.
 
-if nargin < 4 || isempty(model)
-    model = 'inputless';
-end
 
 size_S = size(S);  ndimsx = length(size_S)-2;
 Nt = size_S(end); t = [0:Nt-1]*TR;
@@ -57,21 +52,21 @@
 end
 
 params_all = {'kPL', 'kPB', 'kPA', ...
-	'R1P', 'R1L', 'R1A', 'R1B', ...
-	'S0_L', 'S0_B', 'S0_A', ...
-	'Rinj', 'Tarrival', 'Tbolus'};
+    'R1P', 'R1L', 'R1A', 'R1B', ...
+    'S0_L', 'S0_B', 'S0_A', ...
+    'Rinj', 'Tarrival', 'Tbolus'};
 params_default_est = [0.01, 0.01, 0.01, ...
-	1/30, 1/25, 1/25, 1/15, ...
-	0, 0, 0, ...
-	0.1, 0, 8];
+    1/30, 1/25, 1/25, 1/15, ...
+    0, 0, 0, ...
+    0.1, 0, 8];
 params_default_lb = [-Inf, -Inf, -Inf, ...
-	 1/50, 1/50, 1/50, 1/50, ...
-	 -Inf, -Inf, -Inf, ...
-	  0, -30, 0];
+    1/50, 1/50, 1/50, 1/50, ...
+    -Inf, -Inf, -Inf, ...
+    0, -30, 0];
 params_default_ub = [Inf, Inf, Inf, ...
-	1/10, 1/10, 1/10, 1/5 , ...
-	 Inf, Inf, Inf, ...
-	 Inf 30 Inf];
+    1/10, 1/10, 1/10, 1/5 , ...
+    Inf, Inf, Inf, ...
+    Inf 30 Inf];
 
 if nargin < 5 || isempty(params_fixed)
     params_fixed = struct([]);
@@ -83,11 +78,12 @@
 
 % Supports up to 3 metabolic products (e.g. alanine, lactate, bicarb)
 switch Nmets
-	case 2 % assume pyruvate & lactate
-params_fixed.kPA = 0;  params_fixed.S0_A = 0;
-params_fixed.kPB = 0;  params_fixed.S0_B = 0;
-	case 3 % assume pyruvate & lactate & bicarbonate
-params_fixed.kPA = 0;   params_fixed.S0_A = 0;
+    case 2 % assume pyruvate & lactate
+        params_fixed.kPA = 0;  params_fixed.S0_A = 0;  params_fixed.R1A = 1;
+        params_fixed.kPB = 0;  params_fixed.S0_B = 0;  params_fixed.R1B = 1;
+        
+    case 3 % assume pyruvate & lactate & bicarbonate
+        params_fixed.kPA = 0;   params_fixed.S0_A = 0;  params_fixed.R1A = 1;
 end
 
 
@@ -137,57 +133,62 @@
     disp('==== Computing parameter map ====')
 end
 
-S = reshape(S, [prod(Nx), Nmets, Nt]);  % put all spatial locations in first dimension
+Sreshape = reshape(S, [prod(Nx), Nmets, Nt]);  % put all spatial locations in first dimension
+if Nmets < 4
+    Sreshape = cat(2, Sreshape, zeros([prod(Nx) 4-Nmets, Nt]));  % add zero data for unused metabolites
+    flips = cat(1, flips, ones([4-Nmets, size(flips,2)]));
+end
 
 [Sscale, Mzscale] = flips_scaling_factors(flips, Nt);
 
 params_fit_vec = zeros([prod(Nx),Nparams_to_fit]);  objective_val = zeros([1,prod(Nx)]);
-Sfit = zeros([prod(Nx),Nmets,Nt]); ufit = zeros([prod(Nx),Nt]);
+Sfit = zeros([prod(Nx),Nmets-1,Nt]); ufit = zeros([prod(Nx),Nt]);
 
-for i=1:size(S, 1)
-    if length(Nx) > 1 && plot_flag
+for i=1:size(Sreshape, 1)
+    if prod(Nx) > 1 && plot_flag
         disp([num2str( floor(100*(i-1)/size(S, 1)) ) '% complete'])
     end
     % observed magnetization (Mxy)
-    Mxy = reshape(S(i, :, :), [Nmets, Nt]); % pyr
-
-    if any(y1 ~= 0)
+    Mxy = reshape(Sreshape(i, :, :), [4, Nt]);
+    
+    if any(Mxy(:) ~= 0)
         
         % estimate state magnetization (MZ) based on scaling from RF pulses
-        x1 = y1./Sscale(1, :);
-        x2 = y2./Sscale(2, :);
+        Mz = Mxy./Sscale;
         
         % fit to data
         options = optimoptions(@fminunc,'Display','none','Algorithm','quasi-newton');
         lsq_opts = optimset('Display','none','MaxIter', 500, 'MaxFunEvals', 500);
+        
         switch(fit_method)
             case 'ls'
-                obj = @(var) (x2 - trajectories_inputless(model, var, x1, Mzscale, params_fixed, TR)) .* Sscale(2,:);  % perform least-squares in signal domain
+                obj = @(var) difference_inputless(var, params_fixed, TR, Mzscale, Sscale, Mz, Nmets) ;  % perform least-squares in signal domain
                 [params_fit_vec(i,:),objective_val(i)] = lsqnonlin(obj, params_est_vec, params_lb, params_ub, lsq_opts);
                 
             case 'ml'
-                obj = @(var) negative_log_likelihood_rician_inputless(model, var, x1, x2, Mzscale, params_fixed, TR, noise_level.*(Sscale(2,:).^2));
+                obj = @(var) negative_log_likelihood_rician_inputless(var, params_fixed, TR, Mzscale, Mz, noise_level.*(Sscale).^2, Nmets);
                 [params_fit_vec(i,:), objective_val(i)] = fminunc(obj, params_est_vec, options);
                 
         end
         
-        [Sfit(i,:), ufit(i,:)] = trajectories(params_fit_vec(i,:), x1, Mzscale, params_fixed, TR);
-        Sfit(i,:) = Sfit(i,:)  .* Sscale(2, :);
+        [Mzfit, ufit(i,:)] = trajectories_inputless(params_fit_vec(i,:), params_fixed, TR,  Mzscale, Mz(1,:));
+        
+        Sfit(i,:,:) = Mzfit(2:Nmets,:)  .* Sscale(2:Nmets, :);
         ufit(i,:) = ufit(i,:)  .* Sscale(1, :);
         
         if plot_flag
             % plot of fit for debugging
             figure(99)
             subplot(2,1,1)
-            plot(t, x1, t, x2, t, Sfit(i,:)./ Sscale(2, :),'--', t, ufit(i,:)./ Sscale(1, :), 'k:')
+            plot(t, Mz, t, Mzfit,'--', t, ufit(i,:)./ Sscale(1, :), 'k:')
             xlabel('time (s)')
             ylabel('state magnetization (au)')
             subplot(2,1,2)
-            plot(t, y1, t, y2, t, Sfit(i,:),'--', t, ufit(i,:), 'k:')
+            plot(t, Mxy, t, squeeze(Sfit(i,:,:)),'--', t, ufit(i,:), 'k:')
             xlabel('time (s)')
             ylabel('signal (au)')
-            title(num2str(params_fit_vec(i,1:end-1),4)) % don't display L0_start value
-            legend('pyruvate', 'lactate', 'lactate fit', 'input estimate')
+            title(num2str(params_fit_vec(i,:),2)) % don't display L0_start value
+            %            legend('pyruvate', 'lactate', 'lactate fit', 'input estimate')
             drawnow, pause(0.5)
         end
     end
@@ -210,41 +211,48 @@
     end
     
     
-    Sfit = reshape(Sfit, [Nx, Nt]);
+    Sfit = reshape(Sfit, [Nx, Nmets-1, Nt]);
     ufit = reshape(ufit, [Nx, Nt]);
     objective_val = reshape(objective_val, Nx);
-	if plot_flag
-    	disp('100 % complete')
+    if plot_flag
+        disp('100 % complete')
     end
 end
 
 end
 
-function [ l1 ] = negative_log_likelihood_rician_inputless(params_fit, x1, x2, Mzscale, params_fixed, TR, noise_level)
+function diff_products = difference_inputless(params_fit, params_fixed, TR, Mzscale, Sscale, Mz, Nmets)
+Mzfit = trajectories_inputless(params_fit, params_fixed, TR,  Mzscale, Mz(1,:)) ;
+temp_diff = (Mz - Mzfit) .* Sscale;
+diff_products = temp_diff(2:Nmets,:);
+diff_products = diff_products(:);
+end
+
+function [ l1 ] = negative_log_likelihood_rician_inputless(params_fit, params_fixed, TR, Mzscale, Mz, noise_level, Nmets)
 %FUNCTION NEGATIVE_LOG_LIKELIHOOD_RICIAN Computes log likelihood for
 %    compartmental model with Rician noise
 % noise level is scaled for state magnetization (Mz) domain
 
-N = size(x1,2);
+N = size(Mzscale,2);
 
 % compute trajectory of the model with parameter values
-x2fit = trajectories_inputless(params_fit, x1, Mzscale, params_fixed, TR);
+Mzfit = trajectories_inputless(params_fit, params_fixed, TR,  Mzscale, Mz(1,:)) ;
 
 % compute negative log likelihood
 l1 = 0;
 for t = 1:N
-    for k = 1
+    for k = 2:Nmets
         l1 = l1 - (...
-            log(x2(k, t)) - log(noise_level(t)) ...
-            - (x2(k, t)^2 + x2fit(k, t)^2)/(2*noise_level(t)) ...
-            + x2(k, t)*x2fit(k, t)/noise_level(t) ...
-            + log(besseli(0, x2(k, t)*x2fit(k, t)/noise_level(t), 1))...
+            log(Mz(k, t)) - log(noise_level(k,t)) ...
+            - (Mz(k, t)^2 + Mzfit(k, t)^2)/(2*noise_level(k,t)) ...
+            + Mz(k, t)*Mzfit(k, t)/noise_level(k,t) ...
+            + log(besseli(0, Mz(k, t)*Mzfit(k, t)/noise_level(k,t), 1))...
             );
     end
 end
 end
 
-function [Mz_products, u] = trajectories_inputless( params_fit, Mz_pyr, Mzscale, params_fixed , TR )
+function [Mz_all, u] = trajectories_inputless( params_fit, params_fixed, TR, Mzscale , Mz_pyr )
 % Compute product magnetizations using a uni-directional two-site model
 % Uses substrate magnetization measurements, estimated relaxation and
 % conversion rates
@@ -254,8 +262,8 @@
 u = zeros(1,N);
 
 params_all = {'kPL', 'kPB', 'kPA', ...
-	'R1P', 'R1L', 'R1A', 'R1B', ...
-	'S0_L', 'S0_B', 'S0_A'};
+    'R1P', 'R1L', 'R1A', 'R1B', ...
+    'S0_L', 'S0_B', 'S0_A'};
 
 nfit = 0;
 for n = 1:length(params_all)
@@ -272,15 +280,23 @@
 Mz_all(3,1) = S0_B;
 Mz_all(4,1) = S0_A;
 
+A = [-R1P-kPL-kPB-kPA, 0, 0, 0
+    +kPL, -R1L, 0, 0
+    +kPB, 0, -R1B, 0
+    +kPA, 0, 0, -R1A];
+
 for It=1:N-1
     
     Mz_init = Mz_all(:,It) .* Mzscale(:, It);
     
     % estimate input, assuming this is constant during TR interval
-    u(It) = ( Mz_pyr(It+1) - Mz_init(1)*exp((- R1P - kPL - kPB - kPA)*TR) ) * (R1P + kPL) / (1 - exp((- R1P - kPL - kPB - kPA)*TR));
+    % This calculation could be improved for noise stability?
+    u(It) = ( Mz_pyr(It+1) - Mz_init(1)*exp((- R1P - kPL - kPB - kPA)*TR) ) * (R1P + kPL + kPB + kPA) / (1 - exp((- R1P - kPL - kPB - kPA)*TR));
+    
+    xstar = - inv(A)*[u(It),0,0,0].';
     
     % solve next time point under assumption of constant input during TR
-    x2(It+1) = exp(-R1L*TR)*((L0*R1L*R1P - L0*R1L^2 - kPL*u(It) + L0*R1L*kPL + P0*R1L*kPL)/(R1L*(R1P - R1L + kPL)) + (kPL*u(It)*exp(R1L*TR))/(R1L*(R1P - R1L + kPL))) - exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u(It) + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u(It)*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)));
+    Mz_all(:,It+1) = xstar + expm(A*TR) * (Mz_init - xstar);
     
     
 end
diff --git a/kinetic_modeling/test_fit_pyr_kinetics.m b/kinetic_modeling/test_fit_pyr_kinetics.m
new file mode 100644
index 0000000..e636e28
--- /dev/null
+++ b/kinetic_modeling/test_fit_pyr_kinetics.m
@@ -0,0 +1,196 @@
+% Script for testing fit_kPL kinetic model fitting functions
+
+clear all
+
+% Test values
+Tin = 0; Tacq = 48; TR = 3; N = Tacq/TR;
+R1P = 1/25; R1L = 1/25; R1B = 1/15; R1A = 1/25;
+kPL = 0.05; kPB = 0.03; kPA = 0.02;
+std_noise = 0.005;
+
+% gamma variate input function - most realistic
+t = [1:N]*TR;
+Tarrival = 0;
+A = 4; B = 1*TR;
+input_function = gampdf(t-Tarrival,A,B);  % gamma distribution -continued input
+input_function = input_function/sum(input_function);% normalize for a total magnetization input = 1
+Mz0 = zeros(4,1);
+
+% Test over multiple combinations of flip angle schemes
+k12 = 0.03; % for variable flip angle designs
+flips(1:2,1:N,1) = ones(2,N)*30*pi/180;  % constant, single-band
+flips(1:2,1:N,2) = repmat([20;35]*pi/180,[1 N]);  % constant, multi-band
+flips(1:2,1:N,3) = [vfa_const_amp(N, pi/2, exp(-TR * ( k12))); ... % max lactate SNR variable flip angle
+    vfa_opt_signal(N, exp(-TR * ( R1L)))];
+
+flips = cat(1, flips, repmat( flips(2,:,:), [2 1 1]));
+
+N_flip_schemes = size(flips,3);
+
+t = [0:N-1]*TR + Tin;
+
+figure
+subplot(121) , plot(t, squeeze(flips(1,:,:))*180/pi)
+title('Pyruvate flips')
+subplot(122) , plot(t, squeeze(flips(2,:,:))*180/pi)
+title('Product (Lactate, Alanine, Bicarbonate) flips')
+legend('constant','multiband',  'max lactate SNR vfa'); %, 'vfa', 'T1-effective vfa', 'Saturation Recovery')
+
+% generate simulated data
+noise_S = randn([4 N])*std_noise;  % same noise for all flip schedules
+for Iflips = 1:N_flip_schemes
+    [Mxy(1:4, 1:N, Iflips), Mz] = simulate_Nsite_model(Mz0, [R1P R1L R1B R1A], [kPL 0; kPB 0; kPA 0], flips(:,:,Iflips), TR, input_function);
+    % add noise
+    Sn(1:size(Mxy,1), 1:size(Mxy,2),  Iflips) = Mxy(:,:,Iflips) + noise_S;
+end
+
+% initial parameter guesses
+R1P_est = 1/25; R1L_est = 1/25; R1B_est = 1/15; R1A_est = 1/25;
+kPL_est = .02; kPB_est = .02; kPA_est = .02;
+plot_fits = 1;
+
+%% Test fitting - fit kPX only
+disp('Fitting kPL, kPB, and kPA, with fixed relaxation rates:')
+disp('')
+
+clear params_fixed params_est params_fit params_fitn_complex params_fitn_mag
+params_fixed.R1P = R1P_est; params_fixed.R1L = R1L_est; params_fixed.R1B = R1B_est; params_fixed.R1A = R1A_est;
+params_est.kPL = kPL_est; params_est.kPB = kPB_est; params_est.kPA = kPA_est;
+
+for Iflips = 1:N_flip_schemes
+    % no noise
+    [params_fit(:,Iflips) Sfit(1:3,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(Mxy(:,:,Iflips), TR, flips(:,:,Iflips), params_fixed, params_est, [], plot_fits);
+    
+    % add noise
+    [params_fitn_complex(:,Iflips) Snfit_complex(1:3,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(Sn(:,:,Iflips), TR, flips(:,:,Iflips), params_fixed, params_est, [], plot_fits);
+    
+    % magnitude fitting with noise
+ %   [params_fitn_mag(:,Iflips) Snfit_mag(1:3,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(abs(Sn(:,:,Iflips)), TR, flips(:,:,Iflips),params_fixed, params_est, std_noise, plot_fits);
+end
+
+disp('Input:')
+disp(['KPL  = ' num2str(kPL)])
+disp(['KPB  = ' num2str(kPB)])
+disp(['KPA  = ' num2str(kPA)])
+
+disp('Noiseless fit results:')
+k_fit = struct2array(params_fit);
+k_fit = k_fit([[1:3],[9:11],[17:19]]); 
+disp([['KPL  = ';'KPB  = ';'KPA  = '],num2str(reshape(k_fit,[3 3]),2)])
+disp('Noisy complex fit results:')
+k_fit = struct2array(params_fitn_complex);
+k_fit = k_fit([[1:3],[9:11],[17:19]]); 
+disp([['KPL  = ';'KPB  = ';'KPA  = '],num2str(reshape(k_fit,[3 3]),2)])
+% disp('Noisy magnitude fit results:')
+% k_fit = struct2array(params_fitn_mag);
+% k_fit = k_fit([[1:3],[9:11],[17:19]]); 
+% disp([['KPL  = ';'KPB  = ';'KPA  = '],num2str(reshape(k_fit,[3 3]),2)])
+
+figure
+subplot(221) , plot(t, squeeze(Sn(1,:,:)))
+title('Pyruvate signals')
+subplot(222) , plot(t, squeeze(Sn(2,:,:)))
+hold on, plot(t, squeeze(Snfit_complex(1,:,:)),':')
+%plot(t, squeeze(Snfit_mag(1,:,:)),'--'), hold off
+title('Lactate signals and fits')% (dots=complex fit, dashed=magnitude)')
+legend('constant','multiband', 'multiband variable flip')
+subplot(223) , plot(t, squeeze(Sn(3,:,:)))
+hold on, plot(t, squeeze(Snfit_complex(2,:,:)),':')
+%plot(t, squeeze(Snfit_mag(2,:,:)),'--'), hold off
+title('Bicarb signals and fits')
+subplot(224) , plot(t, squeeze(Sn(4,:,:)))
+hold on, plot(t, squeeze(Snfit_complex(3,:,:)),':')
+%plot(t, squeeze(Snfit_mag(3,:,:)),'--'), hold off
+title('Alanine signals and fits')
+
+return
+%% Test fitting - 3-site
+
+clear params_fixed params_est params_fit params_fitn_complex params_fitn_mag
+clear Sfit Snfit_complex Snfit_mag
+params_fixed.R1P = R1P_est; params_fixed.R1L = R1L_est; params_fixed.R1B = R1B_est; params_fixed.R1A = R1A_est;
+params_est.kPL = kPL_est; params_est.kPB = kPB_est; params_est.kPA = kPA_est;
+
+for Iflips = 1:N_flip_schemes
+    % no noise
+    [params_fit(:,Iflips) Sfit(1:2,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(Mxy(1:3,:,Iflips), TR, flips(1:3,:,Iflips), params_fixed, params_est, [], plot_fits);
+    
+    % add noise
+    [params_fitn_complex(:,Iflips) Snfit_complex(1:2,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(Sn(1:3,:,Iflips), TR, flips(1:3,:,Iflips), params_fixed, params_est, [], plot_fits);
+    
+    % magnitude fitting with noise
+    [params_fitn_mag(:,Iflips) Snfit_mag(1:2,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(abs(Sn(1:3,:,Iflips)), TR, flips(1:3,:,Iflips),params_fixed, params_est, std_noise, plot_fits);
+end
+
+disp('Input:')
+disp(['KPL  = ' num2str(kPL)])
+disp(['KPB  = ' num2str(kPB)])
+
+disp('Noiseless fit results:')
+k_fit = struct2array(params_fit);
+k_fit = k_fit([[1:2],[7:8],[13:14]]); 
+disp([['KPL  = ';'KPB  = '],num2str(reshape(k_fit,[2 3]),2)])
+disp('Noisy complex fit results:')
+k_fit = struct2array(params_fitn_complex);
+k_fit = k_fit([[1:2],[7:8],[13:14]]); 
+disp([['KPL  = ';'KPB  = '],num2str(reshape(k_fit,[2 3]),2)])
+disp('Noisy magnitude fit results:')
+k_fit = struct2array(params_fitn_mag);
+k_fit = k_fit([[1:2],[7:8],[13:14]]); 
+disp([['KPL  = ';'KPB  = '],num2str(reshape(k_fit,[2 3]),2)])
+
+figure
+subplot(221) , plot(t, squeeze(Sn(1,:,:)))
+title('Pyruvate signals')
+subplot(222) , plot(t, squeeze(Sn(2,:,:)))
+hold on, plot(t, squeeze(Snfit_complex(1,:,:)),':')
+plot(t, squeeze(Snfit_mag(1,:,:)),'--'), hold off
+title('Lactate signals and fits (dots=complex fit, dashed=magnitude)')
+legend('constant','multiband', 'multiband variable flip')
+subplot(223) , plot(t, squeeze(Sn(3,:,:)))
+hold on, plot(t, squeeze(Snfit_complex(2,:,:)),':')
+plot(t, squeeze(Snfit_mag(2,:,:)),'--'), hold off
+title('Bicarb signals and fits')
+
+%% Test fitting - 2-site
+
+clear params_fixed params_est params_fit params_fitn_complex params_fitn_mag
+clear Sfit Snfit_complex Snfit_mag
+params_fixed.R1P = R1P_est; params_fixed.R1L = R1L_est; params_fixed.R1B = R1B_est; params_fixed.R1A = R1A_est;
+params_est.kPL = kPL_est; params_est.kPB = kPB_est; params_est.kPA = kPA_est;
+
+for Iflips = 1:N_flip_schemes
+    % no noise
+    [params_fit(:,Iflips) Sfit(1,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(Mxy(1:2,:,Iflips), TR, flips(1:2,:,Iflips), params_fixed, params_est, [], plot_fits);
+    
+    % add noise
+    [params_fitn_complex(:,Iflips) Snfit_complex(1,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(Sn(1:2,:,Iflips), TR, flips(1:2,:,Iflips), params_fixed, params_est, [], plot_fits);
+    
+    % magnitude fitting with noise
+    [params_fitn_mag(:,Iflips) Snfit_mag(1,1:size(Mxy,2),  Iflips)] = fit_pyr_kinetics(abs(Sn(1:2,:,Iflips)), TR, flips(1:2,:,Iflips),params_fixed, params_est, std_noise, plot_fits);
+end
+
+disp('Input:')
+disp(['KPL  = ' num2str(kPL)])
+
+disp('Noiseless fit results:')
+k_fit = struct2array(params_fit);
+k_fit = k_fit([[1],[5],[9]]); 
+disp([['KPL  = '],num2str(reshape(k_fit,[1 3]),2)])
+disp('Noisy complex fit results:')
+k_fit = struct2array(params_fitn_complex);
+k_fit = k_fit([[1],[5],[9]]); 
+disp([['KPL  = '],num2str(reshape(k_fit,[1 3]),2)])
+disp('Noisy magnitude fit results:')
+k_fit = struct2array(params_fitn_mag);
+k_fit = k_fit([[1],[5],[9]]); 
+disp([['KPL  = '],num2str(reshape(k_fit,[1 3]),2)])
+
+figure
+subplot(221) , plot(t, squeeze(Sn(1,:,:)))
+title('Pyruvate signals')
+subplot(222) , plot(t, squeeze(Sn(2,:,:)))
+hold on, plot(t, squeeze(Snfit_complex(1,:,:)),':')
+plot(t, squeeze(Snfit_mag(1,:,:)),'--'), hold off
+title('Lactate signals and fits (dots=complex fit, dashed=magnitude)')
+legend('constant','multiband', 'multiband variable flip')
diff --git a/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m b/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m
index 70cb78f..87d5b64 100644
--- a/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m	
+++ b/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m	
@@ -1,19 +1,18 @@
 % fit slab dynamic MRS data from repeated injections ("shots") in a single rat
-% as an example, does a kPL fit only (ignores alanine, bicarb)
 
 % Data and studies described in
-% Hu, Simon, Peder E Z Larson, Mark Vancriekinge, Andrew M Leach, Ilwoo Park, Christine Leon,  et al. 
-% “Rapid Sequential Injections of Hyperpolarized [1-¹³C]pyruvate in Vivo Using a Sub-Kelvin, Multi-Sample DNP Polarizer.” 
-% Magn Reson Imaging 31, no. 4 (May 2013): 490–96. 
+% Hu, Simon, Peder E Z Larson, Mark Vancriekinge, Andrew M Leach, Ilwoo Park, Christine Leon,  et al.
+% “Rapid Sequential Injections of Hyperpolarized [1-¹³C]pyruvate in Vivo Using a Sub-Kelvin, Multi-Sample DNP Polarizer.”
+% Magn Reson Imaging 31, no. 4 (May 2013): 490–96.
 % https://doi.org/10.1016/j.mri.2012.09.002.
 
 clear all
 plot_flag = 1;
 int_range = 25;
-    R1P_est = 1/30; R1L_est = 1/25;  R1B_est = 1/15; R1A_est = 1/25;
-    kPL_est = 0.01; kPB_est = .01; kPA_est = .01;
+R1P_est = 1/30; R1L_est = 1/25;  R1B_est = 1/15; R1A_est = 1/25;
+kPL_est = 0.01; kPB_est = .01; kPA_est = .01;
 
-for ratnum = 1:4;
+for ratnum = 1:5
     
     switch ratnum
         case 1
@@ -65,10 +64,12 @@
         load(sprintf('rat%d/shot%d',ratnum,n));
         pyr = real(sum(spectra_dynamic(pyr_center+[-int_range:int_range],:),1));
         lac = real(sum(spectra_dynamic(lac_center+[-int_range:int_range],:),1));
-        ala = real(sum(spectra_dynamic(ala+[-int_range:int_range],:),1));
-        bicarb = real(sum(spectra_dynamic(bicarb+[-int_range:int_range],:),1));
-        [params_fit Sfit] = fit_kPL([pyr;lac] , TR, repmat(flip, [2 length(pyr)]), params_fixed, params_est, [], plot_flag);
-        kPL(n) = params_fit.kPL;
+        ala = real(sum(spectra_dynamic(ala_center+[-int_range:int_range],:),1));
+        bicarb = real(sum(spectra_dynamic(bicarb_center+[-int_range:int_range],:),1));
+        
+        % only fit kPL - gives similar results
+        %         [params_fit Sfit] = fit_kPL([pyr;lac] , TR, repmat(flip, [2 length(pyr)]), params_fixed, params_est, [], plot_flag);
+        %         kPL_only(n) = params_fit.kPL;
         
         [params_fit Sfit] = fit_pyr_kinetics([pyr;lac;bicarb;ala] , TR, repmat(flip, [4 length(pyr)]), params_fixed, params_est, [], plot_flag);
         kPL(n) = params_fit.kPL; kPB(n) = params_fit.kPB; kPA(n) = params_fit.kPA;
@@ -76,7 +77,7 @@
     end
     
     disp(['Rat #' int2str(ratnum)]);
-    disp(['   kPL fits: ' num2str(kPL) ' 1/s'])
-    disp(['   kPB fits: ' num2str(kPB) ' 1/s'])
-    disp(['   kPA fits: ' num2str(kPA) ' 1/s'])
+    disp(['   kPL fits: ' num2str(kPL,3) ' 1/s'])
+    disp(['   kPB fits: ' num2str(kPB,3) ' 1/s'])
+    disp(['   kPA fits: ' num2str(kPA,3) ' 1/s'])
 end
\ No newline at end of file
diff --git a/simulations/simulate_Nsite_model.m b/simulations/simulate_Nsite_model.m
new file mode 100644
index 0000000..0daad7a
--- /dev/null
+++ b/simulations/simulate_Nsite_model.m
@@ -0,0 +1,80 @@
+function [Mxy, Mz] = simulate_Nsite_model(Tin, R1, k, flips, TR, input_function)
+% [Mxy, Mz] = simulate_Nsite_model(Mz0, R1, k, flips, TR, [input_function])
+%
+% Simulates the magnetization evolution in a N-site exchange model with
+% single substrate and multiple products for varying flip angles.
+% Allows for up to N=4
+% Applied disretely as flips(:,1), TR, flips(:,2), TR, etc...
+%
+% INPUTS:
+%   Mz0 [1xN] -  model starts from this
+%       magnetization distribution
+%   R1 - relaxation times for N sites (1/s)
+%   k [Nx2] - conversion rates between sites (1/s), first column is forward
+%   rate, second column is reverse rate
+%   flips - flip angles (radians)
+%   TR - repetition time (s)
+%   input_function (optional) - additional input function to add to
+%   pyruvate for all time points.
+%
+% OUTPUTS:
+%   Mxy,Mz - expected magnetization evolution immediately after flips
+%
+% (c) 2013-2014 The Regents of the University of California
+% All Rights Reserved.
+%
+% Author: Peder E. Z. Larson
+
+Nmets = size(flips,1);
+N = size(flips,2);
+
+if isempty(R1)
+    R1 = zeros(1,Nmets);  % infinite relaxation time
+end
+
+if length(R1) == 1
+    R1 = R1*ones(1,Nmets);
+end
+
+if nargin < 6 || isempty(input_function) || all(input_function == 0)
+    use_input_function = 0;
+else
+    use_input_function = 1;
+    Nsim = 100;
+end
+
+switch Nmets
+    case 2
+        A = [-R1(1)-k(1,1) +k(1,2)
+            +k(1,1) -R1(2)-k(1,2)];
+    case 3
+        A = [-R1(1)-k(1,1)-k(2,1) +k(1,2) +k(2,2)
+            +k(1,1) -R1(2)-k(1,2) 0
+            +k(2,1) 0 -R1(3)-k(2,2)];
+    case 4
+        A = [-R1(1)-k(1,1)-k(2,1) +k(1,2) +k(2,2) +k(3,2)
+            +k(1,1) -R1(2)-k(1,2) 0 0
+            +k(2,1) 0 -R1(3)-k(2,2) 0
+            +k(3,1) 0 0 -R1(4)-k(3,2)];
+end
+
+Mz0 = Tin(:);
+
+
+Mxy(1:Nmets,1) = Mz0 .* sin(flips(:,1));
+Mz(1:Nmets,1) = Mz0 .* cos(flips(:,1));
+
+for n = 2:N
+    if use_input_function
+        Mz_m = Mz(:,n-1);
+        % more accurate to spread out input over a number of samples to
+        % avoid unrealistically large signal jumps
+        for ni = 1:Nsim
+            Mz_m = expm(A*TR/Nsim) * (Mz_m + [input_function(n-1)/Nsim;zeros(Nmets-1,1)]);
+        end
+    else
+        Mz_m = expm(A*TR) * (Mz(:,n-1));
+    end
+    Mxy(:,n) =  Mz_m .* sin(flips(:,n));
+    Mz(:,n) = Mz_m .* cos(flips(:,n));
+end
\ No newline at end of file

From b0e1a838830d2e8faddc5a03a9252e377d048113 Mon Sep 17 00:00:00 2001
From: agentmess <peder.larson@ucsf.edu>
Date: Tue, 20 Mar 2018 22:39:39 -0700
Subject: [PATCH 5/6] Minor fix to gammainput

---
 kinetic_modeling/fit_kPL_withgammainput.m     | 20 +++++++++----------
 .../test_fit_kPL_withgammainput.m             |  2 +-
 .../demo_dynamic_mrs_multishot_fit.m          |  0
 3 files changed, 11 insertions(+), 11 deletions(-)
 rename sample_data/{Rat Kidneys Dynamic MRS => }/demo_dynamic_mrs_multishot_fit.m (100%)

diff --git a/kinetic_modeling/fit_kPL_withgammainput.m b/kinetic_modeling/fit_kPL_withgammainput.m
index 493cbd6..880756b 100644
--- a/kinetic_modeling/fit_kPL_withgammainput.m
+++ b/kinetic_modeling/fit_kPL_withgammainput.m
@@ -109,7 +109,7 @@
 Sfit = zeros([prod(Nx),2,Nt]);
 Sfit1 = zeros([prod(Nx),Nt]); Sfit2 = zeros([prod(Nx),Nt]);
 
-parfor i=1:size(S, 1)
+for i=1:size(S, 1)
     if length(Nx) > 1 && plot_flag
         disp([num2str( floor(100*(i-1)/size(S, 1)) ) '% complete'])
     end
@@ -134,7 +134,7 @@
         switch(fit_method)
             case 'ls'
                 % Fit Pyruvate first (Tarrive, Rinj, Tend...) - seems to help a little
-                obj = @(var) trajectory_difference_x1(var, x1, x2, params_fixed, TR, Mzscale);
+                obj = @(var) trajectory_difference_x1(var, x1, x2, params_fixed, TR, Mzscale, Sscale);
                 [params_fit_vec_temp] = lsqnonlin(obj, params_est_vec, params_lb, params_ub, lsq_opts);
                 
                 params_est_vec2 = params_est_vec;
@@ -147,12 +147,12 @@
                 end
                 
                 % Fit all data
-                obj = @(var) trajectory_difference_all(var, x1, x2, params_fixed, TR, Mzscale);  % add Sscale?
+                obj = @(var) trajectory_difference_all(var, x1, x2, params_fixed, TR, Mzscale, Sscale);  
                 [params_fit_vec(i,:),objective_val(i)] = lsqnonlin(obj, params_est_vec2, params_lb, params_ub, lsq_opts);
                 
             case 'ml'
                 obj = @(var) negative_log_likelihood_rician(var, x1, x2, Mzscale, params_fixed, TR, noise_level.*(Sscale(2,:).^2));
-                [params_fit_vec(i,:), objective_val(i)] = fminunc(obj, params_est_vec2, options);
+                [params_fit_vec(i,:), objective_val(i)] = fminunc(obj, params_est_vec, options);
                 
         end
         [x1fit, x2fit] = trajectories_withgammainput(params_fit_vec(i,:), params_fixed, TR, Nt, Mzscale);
@@ -202,19 +202,19 @@
 
 end
 
-function diff_all = trajectory_difference_all(params_fit, x1, x2,  params_fixed, TR, Mzscale)
+function diff_yall = trajectory_difference_all(params_fit, x1, x2,  params_fixed, TR, Mzscale, Sscale)
 [x1fit, x2fit] = trajectories_withgammainput(params_fit, params_fixed, TR, length(x1), Mzscale) ;
-diff_all = [ x1(:)-x1fit(:) ; x2(:)-x2fit(:)];
+diff_yall = [ (x1(:)-x1fit(:)) .* Sscale(1,:).' ; ( x2(:)-x2fit(:) ) .* Sscale(2,:).'];
 end
 
-function diff_x1 = trajectory_difference_x1(params_fit, x1, x2,  params_fixed, TR, Mzscale)
+function diff_y1 = trajectory_difference_x1(params_fit, x1, x2,  params_fixed, TR, Mzscale, Sscale)
 [x1fit, x2fit] = trajectories_withgammainput(params_fit, params_fixed, TR, length(x1), Mzscale) ;
-diff_x1 = [ x1(:)-x1fit(:) ];
+diff_y1 = [ (x1(:)-x1fit(:)) .* Sscale(1,:).'  ];
 end
 
-function diff_x2 = trajectory_difference_x2(params_fit, x1, x2,  params_fixed, TR, Mzscale)
+function diff_y2 = trajectory_difference_x2(params_fit, x1, x2,  params_fixed, TR, Mzscale, Sscale)
 [x1fit, x2fit] = trajectories_withgammainput(params_fit, params_fixed, TR, length(x1), Mzscale) ;
-diff_x2 = [ x2(:)-x2fit(:) ];
+diff_y2 = [ ( x2(:)-x2fit(:) ) .* Sscale(2,:).'];
 end
 
 
diff --git a/kinetic_modeling/test_fit_kPL_withgammainput.m b/kinetic_modeling/test_fit_kPL_withgammainput.m
index 499e060..0c72474 100644
--- a/kinetic_modeling/test_fit_kPL_withgammainput.m
+++ b/kinetic_modeling/test_fit_kPL_withgammainput.m
@@ -54,7 +54,7 @@
 % generate simulated data
 noise_S = randn([2 N])*std_noise;  % same noise for all flip schedules
 for Iflips = 1:N_flip_schemes
-    [Mxy(1:2, 1:N, Iflips), Mz] = simulate_2site_model(Mz0, [R1P R1L], [KPL 0], flips(:,:,Iflips), TR, input_function);
+    [Mxy(1:2, 1:N, Iflips), Mz] = simulate_Nsite_model(Mz0, [R1P R1L], [KPL 0], flips(:,:,Iflips), TR, input_function);
     % add noise
     Sn(1:size(Mxy,1), 1:size(Mxy,2),  Iflips) = Mxy(:,:,Iflips) + noise_S;
 end
diff --git a/sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m b/sample_data/demo_dynamic_mrs_multishot_fit.m
similarity index 100%
rename from sample_data/Rat Kidneys Dynamic MRS/demo_dynamic_mrs_multishot_fit.m
rename to sample_data/demo_dynamic_mrs_multishot_fit.m

From ea99a083671d75004bb94857a46bc2919db400d0 Mon Sep 17 00:00:00 2001
From: agentmess <peder.larson@ucsf.edu>
Date: Tue, 20 Mar 2018 22:43:22 -0700
Subject: [PATCH 6/6] File reorganization

---
 .../demo_dynamic_mrs_multishot_fit.m          |  1 +
 kinetic_modeling/demo_epi_data_fit.m          |  1 +
 kinetic_modeling/trajectories.m               | 79 -------------------
 3 files changed, 2 insertions(+), 79 deletions(-)
 create mode 120000 kinetic_modeling/demo_dynamic_mrs_multishot_fit.m
 create mode 120000 kinetic_modeling/demo_epi_data_fit.m
 delete mode 100644 kinetic_modeling/trajectories.m

diff --git a/kinetic_modeling/demo_dynamic_mrs_multishot_fit.m b/kinetic_modeling/demo_dynamic_mrs_multishot_fit.m
new file mode 120000
index 0000000..3d8845b
--- /dev/null
+++ b/kinetic_modeling/demo_dynamic_mrs_multishot_fit.m
@@ -0,0 +1 @@
+../sample_data/demo_dynamic_mrs_multishot_fit.m
\ No newline at end of file
diff --git a/kinetic_modeling/demo_epi_data_fit.m b/kinetic_modeling/demo_epi_data_fit.m
new file mode 120000
index 0000000..87f1f55
--- /dev/null
+++ b/kinetic_modeling/demo_epi_data_fit.m
@@ -0,0 +1 @@
+../sample_data/demo_epi_data_fit.m
\ No newline at end of file
diff --git a/kinetic_modeling/trajectories.m b/kinetic_modeling/trajectories.m
deleted file mode 100644
index 37f3dd1..0000000
--- a/kinetic_modeling/trajectories.m
+++ /dev/null
@@ -1,79 +0,0 @@
-function [Mz u] = trajectories( model, params_fit, params_fixed , TR, Mzscale )
-
-% all using MZ component so can account for variable flip angles
-Nmets = size(Mzscale,1); N = size(Mzscale,2);
-Mz = zeros(Nmets, N);
-
-
-params_all = {'kPL', 'kPB', 'kPA', ...
-	'R1P', 'R1L', 'R1A', 'R1B', ...
-	'S0_L', 'S0_B', 'S0_A', ...
-	'Rinj', 'Tarrival', 'Tbolus'};
-
-nfit = 0;
-for n = 1:length(params_all)
-    if isfield(params_fixed, params_all(n))
-        eval([params_all{n} '= params_fixed.(params_all{n});']);
-    else
-        nfit = nfit+1;
-        eval([params_all{n} '= params_fit(nfit);']);
-    end
-end
-
-% gamma input parameters
-A = Tbolus;
-B = Tbolus/4; % assumption that leads to reasonable looking shape
-
-% evolution matrix
-switch Nmets
-    case 2
-        Amat = [-R1P - kPL, 0; ...
-            kPL, -R1L];
-    case 3
-        Amat = [-R1P - kPL - kPB, 0, 0; ...
-            kPL, -R1L, 0; ...
-            kPB, 0, -R1B];
-    case 4
-        Amat = [-R1P - kPL - kPB - kPA, 0, 0, 0; ...
-            kPL, -R1L, 0, 0; ...
-            kPB, 0, -R1B, 0; ...
-            kPA, 0, 0, -R1A];
-end
-
-for It=1:N
-    
-    if strcmp(model, 'inputless')
-        Mz_init = Mz(:,It-1) .* Mzscale(:, It-1);
-        % FIX:
-         u(It) = ( x1(t+1) - P0*exp((- R1P - kPL)*TR) ) * (R1P + kPL) / (1 - exp((- R1P - kPL)*TR));
-    else
-        
-    t = (It-1)*TR;  % solving for Mz at time t (accounting for previous TR interval)
-    
-    if It == 1
-    	% initial Mz in TR after RF pulse
-        Mz_init = zeros(NMets,1);
-
-        if Tarrival < 0
-            % account for longer period of bolus prior to imaging
-            t_preimaging = [floor(Tarrival/TR):0]*TR; % t = 0
-           u(It) = sum( gampdf(t_preimaging-Tarrival,A,B)*Rinj );
-        else
-           u(It) = gampdf(t-Tarrival,A,B)*Rinj;
-        end
-        
-    else
-    	% initial Mz in TR after RF pulse
-    	Mz_init = Mz(:,It-1) .* Mzscale(:, It-1);
-
-        u(It) = gampdf(t-Tarrival,A,B)*Rinj;
-    end
-    end
-        
-    % solve next time point under assumption of constant input during TR
-	Mz(:,It) = expm(Amat*TR)*Mz_init + u; % WRONG WRONG
-	    
-end
-
-end
-