Minot updates and new fitting that ignores input during TR

kinetic_modeling/fit_kPL.m

@@ -123,7 +123,7 @@
     % observed magnetization (Mxy)
     y1 = reshape(S(i, 1, :), [1, Nt]); % pyr
     y2 = reshape(S(i, 2, :), [1, Nt]); % lac
-    if any(y1 ~= 0)
+    if any(y1(:) ~= 0)
         % % plot of observed data for debugging
         % figure(1)
         % plot(t, y1, t, y2)
@@ -132,8 +132,8 @@
         % legend('pyruvate', 'lactate')
 
         % estimate state magnetization (MZ) based on scaling from RF pulses
-        x1 = y1./Sscale(1, :);
-        x2 = y2./Sscale(2, :);
+        x1 = y1./(Sscale(1, :)+eps);  % add eps to eliminate divide by zero errors
+        x2 = y2./(Sscale(2, :)+eps);
 
         % fit to data
         options = optimoptions(@fminunc,'Display','none','Algorithm','quasi-newton');

kinetic_modeling/fit_kPL_noinput.m

@@ -0,0 +1,265 @@
+function [params_fit, Sfit, ufit, objective_val] = fit_kPL_noinput(S, TR, flips, params_fixed, params_est, noise_level, plot_flag)
+% fit_kPL - Simple kinetic model fitting 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 (radians) - 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 curve for lactate
+%   ufit - derived input function (unitless)
+%   objective_val - measure of fit error
+%
+% EXAMPLES - see test_fit_kPL_fcn.m
+%
+% Authors: John Maidens,  Peder E. Z. Larson
+%
+% (c)2015-2017 The Regents of the University of California. All Rights
+% Reserved.
+
+params_all = {'kPL', 'R1L', 'R1P', 'L0_start'};
+params_default_est = [0.01, 1/25, 1/25, 0];
+params_default_lb = [-Inf, 1/60, 1/60, 0];
+params_default_ub = [Inf, 1/10, 1/10, Inf];
+
+if nargin < 4 || isempty(params_fixed)
+    params_fixed = struct([]);
+end
+
+if nargin < 5 || 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 
+    noise_level = [];
+end
+
+if 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]);
+
+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, :)+eps);  % add eps to eliminate divide by zero errors
+        x2 = y2./(Sscale(2, :)+eps);
+
+        % 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_noinput(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,:) = trajectories_frompyr_noinput(params_fit_vec(i,:), x1, Mzscale, params_fixed, TR);
+        Sfit(i,:) = Sfit(i,:)  .* Sscale(2, :);
+
+        if plot_flag
+            % plot of fit for debugging
+            figure(99)
+            subplot(2,1,1)
+            plot(t, x1, t, x2, t, Sfit(i,:)./ Sscale(2, :),'--')
+            xlabel('time (s)')
+            ylabel('state magnetization (au)')
+            subplot(2,1,2)
+            plot(t, y1, t, y2, t, Sfit(i,:),'--')
+            xlabel('time (s)')
+            ylabel('signal (au)')
+            title(num2str(params_fit_vec(i,:),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_noinput(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
+
+function x2 = trajectories_frompyr_noinput( params_fit, x1, Mzscale, params_fixed , TR )
+% Compute product magnetization (e.g. lactate) using a uni-directional two-site model
+% Uses substrate magnetization measurements, estimated relaxation and
+% conversion rates
+% x1 and x2 are pyruvate and lactate, respectively, longitudinal magnetization (MZ) component estimates in order to account for variable flip angles
+
+N = length(x1);
+
+x2 = zeros(1, N);
+
+params_all = {'kPL', 'R1L', 'R1P', 'L0_start'};
+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
+
+x2(1) = L0_start;
+
+for t=1:N-1
+
+    P0 = x1(t)*Mzscale(1, t);
+    L0 = x2(t)*Mzscale(2, t);
+
+
+    % solve next time point under assumption of constant input during TR
+     x2(t+1) = exp(-R1L*TR)*((L0*R1L*R1P - L0*R1L^2 + L0*R1L*kPL + P0*R1L*kPL)/(R1L*(R1P - R1L + kPL)) ) - exp(-TR*(R1P + kPL))*((kPL*(P0*R1P + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) );
+
+    % % solution for next source (pyruvate) time point if needed:
+    % x1(t+1) = (exp(-TR*(R1P + kPL))*((kPL*(P0*R1P - u(t) + P0*kPL))/((R1P + kPL)*(R1P - R1L + kPL)) + (kPL*u(t)*exp(R1P*TR + kPL*TR))/((R1P + kPL)*(R1P - R1L + kPL)))*(R1P - R1L + kPL))/kPL;
+
+    % % Solution without an input function:
+    %  x2(t+1) = (-(kPL*exp((- R1P - kPL)*TR) - kPL*exp(-R1L*TR))/(R1P - R1L + kPL))*Mzscale(1, t)*x1(t) ...
+    %      +  exp(-R1L*TR)*Mzscale(2, t)*x2(t);
+
+end
+
+end
+

simulations/HP_montecarlo_evaluation.m

@@ -365,7 +365,7 @@
 function [kPL_fit, AUC_fit] = fitting_simulation(fit_fcn, Mxy, TR, flips, NMC, std_noise, params_fixed, params_est);
 
 kPL_fit = zeros(1,NMC); AUC_fit = zeros(1,NMC);
-parfor n = 1:NMC
+for n = 1:NMC
     Sn = Mxy + randn(size(Mxy))*std_noise;
 
     params_fit = fit_fcn(Sn, TR, flips, params_fixed, params_est, [], 0);