diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..3e3461ae --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.asv diff --git a/SIRsimulator.m b/SIRsimulator.m index ccb271bf..3cb26ef9 100644 --- a/SIRsimulator.m +++ b/SIRsimulator.m @@ -19,6 +19,7 @@ % prob_stay: the probability of staying in the same region per unit time (0.5) % trans_rate: a scalar value, controlling the baseline infectivity + %% output parameters % Rnor_all: A N_regions * T_total matrix, recording the number of normal % alpha-syn in regions @@ -28,9 +29,11 @@ % could be memory-consuming) % Pmis_all: a N_regions * N_regions * T_total matrix, recording the number of misfolded alpha-syn in paths % could be memory-consuming) -% Rnor0: a N_Regions * 1 vector, the population of normal agents in regions before pathogenic spreading -% Pnor0: a N_Regions * 1 vecotr, the population of normal agents in edges before pathogenic spreading +% Rnor0: a N_Regions * 1 vector, the population of normal agents in regions before pathogenic spreading +% Pnor0: a N_Regions * 1 vecotr, the population of normal agents in edges before pathogenic spreading + +%% % make sure the diag is zero sconnDen(eye(N_regions)==1) = 0; sconnLen(eye(N_regions)==1) = 0; @@ -42,9 +45,10 @@ % multinomial distribution % element (i,j) is the probability of moving from region i to edge (i,j) weights = weights ./ repmat(sum(weights, 2), 1, N_regions); +weights(eye(N_regions, 'logical')) = 0; % convert gene expression scores to probabilities -clearance_rate = normcdf(zscore(GBA)); +clearance_rate = normcdf(zscore(GBA)); synthesis_rate = normcdf(zscore(SNCA)); % store the number of normal/misfoled alpha-syn at each time step @@ -55,94 +59,85 @@ [Rnor, Rmis] = deal(zeros(N_regions, 1)); % number of normal/misfolded alpha-syn in regions [Pnor, Pmis] = deal(zeros(N_regions)); % number of normal/misfolded alpha-syn in paths -%% normal alpha-syn growth +% simplification of variables +alphaTerm = (synthesis_rate .* syn_control) .* dt; +betaTerm = exp(-clearance_rate.*dt); +sTerm = 1 ./ sconnLen .* dt .* v; sTerm(isinf(sTerm)) = 0; +wTerm = weights .* dt; +gamma0 = 1 .* trans_rate ./ ROIsize .* dt ; % the probability of getting misfolded + + +%% normal alpha-syn growth % fill the network with normal proteins iter_max = 1000000000; -display('normal alpha synuclein growth') -for t = 1:iter_max +disp('normal alpha synuclein growth'); +for t = 1:iter_max %%% moving process % regions towards paths % movDrt stores the number of proteins towards each region. i.e. % element in kth row lth col denotes the number of proteins in region k % moving towards l - movDrt = repmat(Rnor, 1, N_regions) .* weights; - movDrt = movDrt .* dt; - movDrt(eye(N_regions)==1) = 0; - + movDrt = Rnor .* wTerm; % IMPLICIT EXPANSION + % paths towards regions % update moving - movOut = Pnor .* v ./ sconnLen ; % longer path & smaller v = lower probability of moving out of paths - movOut(sconnLen == 0) = 0; - - Pnor = Pnor - movOut .* dt + movDrt; - Pnor(eye(N_regions)==1) = 0; - + movOut = Pnor .* sTerm; % longer path & smaller v = lower probability of moving out of paths + + Pnor = Pnor - movOut + movDrt; Rtmp = Rnor; - Rnor = Rnor + sum(movOut, 1)' .* dt - sum(movDrt, 2); - + Rnor = Rnor + sum(movOut, 1)' - sum(movDrt, 2); + %%% growth process - Rnor = Rnor - Rnor.* (1-exp(-clearance_rate.*dt)) + (synthesis_rate .* syn_control).*dt ; - - if abs(Rnor - Rtmp) < (1e-7* Rtmp) - break; - end + Rnor = Rnor.*betaTerm + alphaTerm; + + if abs(Rnor - Rtmp) < (1e-7 * Rtmp); break; end end -Pnor0 = Pnor; -Rnor0 = Rnor; %% misfolded protein spreading process +Pnor0 = Pnor; +Rnor0 = Rnor; % inject misfolded alpha-syn Rmis(seed) = init_number; -display('misfolded alpha synuclein spreading') +disp('misfolded alpha synuclein spreading'); for t = 1:T_total %%% moving process % normal proteins: region -->> paths - movDrt_nor = repmat(Rnor, 1, N_regions) .* weights .* dt; - movDrt_nor(eye(N_regions) == 1) = 0; - + movDrt_nor = Rnor .* wTerm; % IMPLICIT EXPANSION + % normal proteins: paths -->> regions - movOut_nor = Pnor .* v ./ sconnLen; - movOut_nor(sconnLen == 0) = 0; - - + movOut_nor = Pnor .* sTerm; + % misfolded proteins: region -->> paths - movDrt_mis = repmat(Rmis, 1, N_regions) .* weights .* dt; - movDrt_mis(eye(N_regions) == 1) = 0; - + movDrt_mis = Rmis .* wTerm; % IMPLICIT EXPANSION + % misfolded proteins: paths -->> regions - movOut_mis = Pmis .* v ./ sconnLen; - movOut_mis(sconnLen == 0) = 0; - + movOut_mis = Pmis .* sTerm; + % update regions and paths - Pnor = Pnor - movOut_nor.* dt + movDrt_nor; - Pnor(eye(N_regions)==1) = 0; - Rnor = Rnor + sum(movOut_nor, 1)'.* dt - sum(movDrt_nor, 2); - - Pmis = Pmis - movOut_mis.*dt + movDrt_mis; - Pmis(eye(N_regions)==1) = 0; - Rmis = Rmis + sum(movOut_mis, 1)'.*dt - sum(movDrt_mis, 2); - - Rnor_cleared = Rnor .* (1-exp(-clearance_rate.* dt)) ; - Rmis_cleared = Rmis .* (1-exp(-clearance_rate.* dt)) ; - % the probability of getting misfolded - gamma0 = 1 .* trans_rate ./ROIsize; - misProb = 1 - exp( -Rmis .* gamma0 .* dt ) ; % trans_rate: default + Pnor = Pnor - movOut_nor + movDrt_nor; + Rnor = Rnor + sum(movOut_nor, 1)' - sum(movDrt_nor, 2); + + Pmis = Pmis - movOut_mis + movDrt_mis; + Rmis = Rmis + sum(movOut_mis, 1)' - sum(movDrt_mis, 2); + + misProb = 1 - exp( -Rmis .* gamma0 ) ; % trans_rate: default % number of newly infected N_misfolded = Rnor .* exp(-clearance_rate) .* misProb ; + % update - Rnor = Rnor - Rnor_cleared - N_misfolded + (synthesis_rate .* syn_control).*dt; - Rmis = Rmis - Rmis_cleared + N_misfolded; + Rnor = Rnor .* betaTerm + alphaTerm - N_misfolded; + Rmis = Rmis .* betaTerm + N_misfolded; Rnor_all(:, t) = Rnor ; Rmis_all(:, t) = Rmis ; - + % uncomment the following lines if you want outputs of alpha-syn in % paths %Pnor_ave(:, :, t) = Pnor; %Pmis_ave(:, :, t) = Pmis; - + end end diff --git a/SIRsimulator2.m b/SIRsimulator2.m new file mode 100644 index 00000000..89b4d3b5 --- /dev/null +++ b/SIRsimulator2.m @@ -0,0 +1,176 @@ +% SIRsimulator.m +function [Rnor_all, Rmis_all, Rnor0, Pnor0, Pnor_all, Pmis_all] = SIRsimulator2(N_regions, v, dt, T_total, GBA, SNCA, sconnLen, sconnDen, ROIsize, seed, syn_control, init_number, prob_stay, trans_rate) +% A function to simulate the spread of misfolded alpha-syn + +%% input parameters (inside parenthesis are values used in the paper) +% N_regions: number of regions (42) +% v: speed (1) +% dt: time step (0.01) +% T_total: total time steps (10000) +% GBA: GBA gene expression (zscore, N_regions * 1 vector) (empirical GBA expression) +% SNCA: SNCA gene expression after normalization (zscore, N_regions * 1 vector) (empirical SNCA expression) +% sconnLen: structural connectivity matrix (length) (estimated from HCP data) +% sconnDen: structural connectivity matrix (strength) (estimated from HCP data) +% ROIsize: region sizes (voxel counts) +% seed: seed region of misfolded alpha-syn injection (choose as you like? (^?^)= here substantia nigra) +% syn_control: a parameter to control the number of voxels in which +% alpha-syn may get synthesized (region size, i.e., ROIsize) +% init_number: number of injected misfolded alpha-syn (1) +% prob_stay: the probability of staying in the same region per unit time (0.5) +% trans_rate: a scalar value, controlling the baseline infectivity + + +%% output parameters +% Rnor_all: A N_regions * T_total matrix, recording the number of normal +% alpha-syn in regions +% Rmis_all: A N_regions * T_total matrix, recording the number of +% misfolded alph-syn in regions +% Pnor_all: a N_regions * N_regions * T_total matrix, recording the number of normal alpha-syn in paths +% could be memory-consuming) +% Pmis_all: a N_regions * N_regions * T_total matrix, recording the number of misfolded alpha-syn in paths +% could be memory-consuming) +% Rnor0: a N_Regions * 1 vector, the population of normal agents in regions before pathogenic spreading +% Pnor0: a N_Regions * 1 vecotr, the population of normal agents in edges before pathogenic spreading + + +%% +% make sure the diag is zero +sconnDen(eye(N_regions)==1) = 0; +sconnLen(eye(N_regions)==1) = 0; + +% set the mobility pattern +weights = sconnDen; +weights = (1 - prob_stay) .* weights + prob_stay .* diag(sum(weights, 2)) ; + +% multinomial distribution +% element (i,j) is the probability of moving from region i to edge (i,j) +weights = weights ./ repmat(sum(weights, 2), 1, N_regions); +weights(eye(N_regions, 'logical')) = 0; + +% convert gene expression scores to probabilities +clearance_rate = normcdf(zscore(GBA)); +synthesis_rate = normcdf(zscore(SNCA)); + +% store the number of normal/misfoled alpha-syn at each time step +[Rnor_all, Rmis_all] = deal( zeros([N_regions, T_total]) ); +[Pnor_all, Pmis_all] = deal( zeros([N_regions, N_regions, T_total]) ); + +% Rnor, Rmis, Pnor, Pmis store results of single simulation at each time +[Rnor, Rmis] = deal(zeros(N_regions, 1)); % number of normal/misfolded alpha-syn in regions +[Pnor, Pmis] = deal(zeros(N_regions)); % number of normal/misfolded alpha-syn in paths + +% simplification of variables +alphaTerm = (synthesis_rate .* syn_control) .* dt; +betaTerm = exp(-clearance_rate.*dt); +sTerm = 1 ./ sconnLen .* dt .* v; sTerm(isinf(sTerm)) = 0; +wTerm = weights .* dt; +gamma0 = 1 .* trans_rate ./ ROIsize .* dt ; % the probability of getting misfolded + + +%% normal alpha-syn growth +% fill the network with normal proteins +iter_max = 1000000000; +disp('normal alpha synuclein growth'); +for t = 1:iter_max + %%% moving process + % regions towards paths + % movDrt stores the number of proteins towards each region. i.e. + % element in kth row lth col denotes the number of proteins in region k + % moving towards l + movDrt = Rnor .* wTerm; % IMPLICIT EXPANSION + + % paths towards regions + % update moving + movOut = Pnor .* sTerm; % longer path & smaller v = lower probability of moving out of paths + + Pnor = Pnor - movOut + movDrt; + Rtmp = Rnor; + Rnor = Rnor + sum(movOut, 1)' - sum(movDrt, 2); + + %%% growth process + Rnor = Rnor.*betaTerm + alphaTerm; % corresponds to Rtest2 +% Rnor = Rnor.*betaTerm + alphaTerm + sum(movOut, 1)' - sum(movDrt, 2); % corresponds to Rtest ; + + if abs(Rnor - Rtmp) < (1e-7 * Rtmp); break; end +end + + +%% +t = tic; +C = eye(N_regions) - diag(betaTerm) + diag(sum(wTerm,2)); Rtest = (C-wTerm)\alphaTerm; +C = diag(1./betaTerm) - eye(N_regions) + diag(sum(wTerm,2)); Rtest2 = (C-wTerm)\(alphaTerm./betaTerm); +Ptest = (Rtest .* wTerm) ./ sTerm; Ptest(sTerm == 0) = 0; +fprintf('Time to estimate closed form for phase 1: %f\n', toc(t)); + +figure; subplot(2, 2, 1); scatter(Rnor, Rtest, [], 1:N_regions, 'filled'); lsline; +xlabel('Numerical solution, R_{nor}'); ylabel('Closed form prediction, R_{test}'); +title('Closed form phase 1 - MG first try'); + +subplot(2, 2, 2); scatter(Rtest, Rtest2, [], 1:N_regions, 'filled'); lsline; +title('Changing order of updating'); +xlabel('R_{test} v1'); ylabel('R_{test} v2'); + +subplot(2, 2, 3); +scatter(Rtest, (Rtest + sum(Ptest .* sTerm,1)' - sum(Rtest.*wTerm,2)).*betaTerm + alphaTerm, [], 1:N_regions, 'filled'); +title("Results of incrementing approximation"); +xlabel("R_{test}^{t}"); ylabel("R_{test}^{t+1}"); + +subplot(2, 2, 4); +scatter(Rnor, (Rnor + sum(Pnor .* sTerm,1)' - sum(Rnor.*wTerm,2)).*betaTerm + alphaTerm, [], 1:N_regions, 'filled'); +xlabel("R_{nor}^{t}"); ylabel("R_{nor}^{t+1}"); + + +%% misfolded protein spreading process + +% % % % % % % % % % % % % % % % +% Rnor = Rtest; Pnor = Ptest; % to use values from closed form calculation +% % % % % % % % % % % % % % % % +Pnor0 = Pnor; +Rnor0 = Rnor; + +% inject misfolded alpha-syn +Rmis(seed) = init_number; +disp('misfolded alpha synuclein spreading'); +for t = 1:T_total + %%% moving process + % normal proteins: region -->> paths + movDrt_nor = Rnor .* wTerm; % IMPLICIT EXPANSION + + % normal proteins: paths -->> regions + movOut_nor = Pnor .* sTerm; + + % misfolded proteins: region -->> paths + movDrt_mis = Rmis .* wTerm; % IMPLICIT EXPANSION + + % misfolded proteins: paths -->> regions + movOut_mis = Pmis .* sTerm; + + % update regions and paths + Pnor = Pnor - movOut_nor + movDrt_nor; + Rnor = Rnor + sum(movOut_nor, 1)' - sum(movDrt_nor, 2); + + Pmis = Pmis - movOut_mis + movDrt_mis; + Rmis = Rmis + sum(movOut_mis, 1)' - sum(movDrt_mis, 2); + + misProb = 1 - exp( -Rmis .* gamma0 ) ; % trans_rate: default + % number of newly infected + N_misfolded = Rnor .* exp(-clearance_rate) .* misProb ; + + % update + Rnor = Rnor .* betaTerm + alphaTerm - N_misfolded; + Rmis = Rmis .* betaTerm + N_misfolded; + + Rnor_all(:, t) = Rnor ; + Rmis_all(:, t) = Rmis ; + + % uncomment the following lines if you want outputs of alpha-syn in + % paths + %Pnor_ave(:, :, t) = Pnor; + %Pmis_ave(:, :, t) = Pmis; + +end +end + + + + diff --git a/main.m b/main.m index ad6f6f7a..cd592c04 100644 --- a/main.m +++ b/main.m @@ -4,7 +4,7 @@ % load gene expressions, real atrophy, ROIsize, functional connectivity... load('data/42regions/workspace.mat'); % load structural connectivity -load('data/42regions/sc35.mat'); +load('data/42regions/sc30.mat'); N_regions = 42; v = 1; @@ -17,8 +17,43 @@ seed = N_regions; % load your GBA, SNCA, sconnDen, sconnLen, ROISize .... +allclose = @(x, y, tol) all(abs(x-y)>> -[Rnor_all, Rmis_all, Rnor0] = SIRsimulator(N_regions, v, dt, T_total, GBA, SNCA, sconnLen, sconnDen, ROIsize, seed, syn_control, init_number, prob_stay, trans_rate); + +f1 = @() SIRsimulator(N_regions, v, dt, T_total, GBA, SNCA, sconnLen, sconnDen, ROIsize, seed, syn_control, init_number, prob_stay, trans_rate); +f2 = @() SIRsimulator2(N_regions, v, dt, T_total, GBA, SNCA, sconnLen, sconnDen, ROIsize, seed, syn_control, init_number, prob_stay, trans_rate); + +clc +tic +[Rnor_all, Rmis_all, Rnor0] = f1(); +toc + +tic +[Rnor_all_2, Rmis_all_2, Rnor0_2] = f2(); +toc + +% timeit(f1), timeit(f2) + +allclose(Rnor_all, Rnor_all_2, 1e-10) +allclose(Rmis_all, Rmis_all_2, 1e-10) +allclose(Rnor0, Rnor0_2, 1e-10) + +%% +figure; +subplot(2, 4, 1); imagesc(Rnor_all'); title('R_{nor}'); +subplot(2, 4, 2); imagesc(Rnor_all_2'); title('R_{nor}^{test}'); +subplot(2, 4, 3); imagesc(Rmis_all'); title('R_{mis}') +subplot(2, 4, 4); imagesc(Rmis_all_2'); title('R_{mis}^{test}'); + +subplot(2, 4, [5 6]); imagesc(abs(Rnor_all' - Rnor_all_2')./Rnor_all'<0.1); +title("Where is R_{nor}^{test} within 10% of expected?") + +subplot(2, 4, [7 8]); imagesc(abs(Rmis_all' - Rmis_all_2')./Rmis_all'<0.1); +title("Where does R_{mis}^{test} within 10% of expected?") + +%% ratio = Rmis_all ./(Rnor_all + Rmis_all) ; ratio(Rmis_all<1) = 0; % remove possible NaNs... diff --git a/mockData.mlx b/mockData.mlx new file mode 100644 index 00000000..3c0b18c0 Binary files /dev/null and b/mockData.mlx differ