Skip to content

Speed up code #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
magnesium2400 wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Speed up code #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6b5e681
testing faster simulator
magnesium2400 Aug 27, 2023
f5d7c92
cleaned up and sped up code. no changes to function (when tested on o…
magnesium2400 Aug 27, 2023
dfcb003
testing closed form solutions
magnesium2400 Aug 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
*.asv
107 changes: 51 additions & 56 deletions SIRsimulator.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
Expand All @@ -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
Expand All @@ -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

Expand Down
176 changes: 176 additions & 0 deletions SIRsimulator2.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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




Loading