From 34540323f37831390db152f14a4e5a3a63eb71c3 Mon Sep 17 00:00:00 2001 From: agentmess Date: Fri, 23 Mar 2018 09:55:05 -0700 Subject: [PATCH] Minot updates and new fitting that ignores input during TR --- kinetic_modeling/fit_kPL.m | 6 +- kinetic_modeling/fit_kPL_noinput.m | 265 +++++++++++++++++++++++++ simulations/HP_montecarlo_evaluation.m | 2 +- 3 files changed, 269 insertions(+), 4 deletions(-) create mode 100644 kinetic_modeling/fit_kPL_noinput.m diff --git a/kinetic_modeling/fit_kPL.m b/kinetic_modeling/fit_kPL.m index 9fcc4ab..7a688dc 100644 --- a/kinetic_modeling/fit_kPL.m +++ b/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'); diff --git a/kinetic_modeling/fit_kPL_noinput.m b/kinetic_modeling/fit_kPL_noinput.m new file mode 100644 index 0000000..de7e50b --- /dev/null +++ b/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 + diff --git a/simulations/HP_montecarlo_evaluation.m b/simulations/HP_montecarlo_evaluation.m index 6aa6ebe..3457a3b 100644 --- a/simulations/HP_montecarlo_evaluation.m +++ b/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);