Alterations in brain spatiotemporal autocorrelation and functional connectivity in lifelong premature ejaculation: A resting-state fMRI study.

1. # Part of the spatiotemporal package for python
2. # Copyright 2022 Max Shinn <[email hidden]>
3. # Available under the MIT license
4. import numpy as np
5. import scipy.signal
6. import scipy.optimize
7. from .tools import spatial_exponential_floor
9. def spatiotemporal_model_timeseries(distance_matrix, sa_lambda, sa_inf, ta_delta1s, num_timepoints, sample_rate, highpass_freq, seed=None):
10. """Simulate the spatiotemporal model from [Shinn et al (2022)](https://www.biorxiv.org/content/10.1101/2021.06.01.446561v1)
12. Args:
13. distance_matrix (NxN numpy array): the NxN distance matrix, representing the spatial distance
14. between location of each of the timeseries. This should usually be the
15. output of the `distance_matrix_euclidean` function.
16. sa_lambda (float): the SA-λ parameter
17. sa_inf (float): the SA-∞ parameter
18. ta_delta1s (length-N list of floats): a list of TA-Δ₁ values, of length N, for generating the timeseries
19. num_timepoints (int): length of timeseries to generate
20. sample_rate (float): the spacing between timepoints (e.g. TR in fMRI)
21. highpass_freq (float): if non-zero, apply a highpass filter above the
22. given frequency. A good default is 0.01 for fMRI timeseries.
23. seed (int, optional): the random seed. If not specified, it will use
24. the current state of the numpy random number generator.
26. Returns:
27. NxT numpy array: For each of the N nodes, a timeseries of length T=`num_timepoints` according to the spatiotemporal model
28. """
29. assert num_timepoints % 2 == 0, "Must be even timeseries length"
30. # Filtered brown noise spectrum
31. spectrum = make_spectrum(tslen=num_timepoints, sample_rate=sample_rate, alpha=2, highpass_freq=highpass_freq)
32. spectra = np.asarray([spectrum]*len(ta_delta1s))
33. # Spatial autocorrelation matrix
34. corr = spatial_exponential_floor(distance_matrix, sa_lambda, sa_inf)
35. # Create spatially embedded timeseries with the given spectra
36. tss = spatial_temporal_timeseries(corr, spectra, seed=seed)
37. # Compute the standard deviation of nosie we need to add to get the desired TA-delta1
38. noises = [how_much_noise(spectrum, max(.001, ta_delta1)) for ta_delta1 in ta_delta1s]
39. # Add noise to the timeseries
40. rng = np.random.RandomState(seed)
41. tss += rng.randn(tss.shape[0], tss.shape[1]) * np.asarray(noises).reshape(-1,1)
42. return tss
44. def spatiotemporal_noiseless_model_timeseries(distance_matrix, sa_lambda, sa_inf, ta_delta1s, num_timepoints, sample_rate, highpass_freq, seed=0):
45. """Simulate the noiseless spatiotemporal model from [Shinn et al (2022)](https://www.biorxiv.org/content/10.1101/2021.06.01.446561v1)
47. Args:
48. distance_matrix (NxN numpy array): the NxN distance matrix, representing the spatial distance
49. between location of each of the timeseries. This should usually be the
50. output of the `distance_matrix_euclidean` function.
51. sa_lambda (float): the SA-λ parameter
52. sa_inf (float): the SA-∞ parameter
53. ta_delta1s (length-N list of floats): a list of TA-Δ₁ values, of length N, for generating the timeseries
54. num_timepoints (int): length of timeseries to generate
55. sample_rate (float): the spacing between timepoints (e.g. TR in fMRI)
56. highpass_freq (float): if non-zero, apply a highpass filter above the
57. given frequency. A good default is 0.01 for fMRI timeseries.
58. seed (int, optional): the random seed. If not specified, it will use
59. the current state of the numpy random number generator.
61. Returns:
62. NxT numpy array: For each of the N nodes, a timeseries of length T=`num_timepoints` according to the noiseless spatiotemporal model
63. """
64. assert num_timepoints % 2 == 0, "Must be even timeseries length"
65. # Determine the pink noise exponent alpha from the TA-delta1
66. alphas = [ta_to_alpha_fast(sample_rate=sample_rate, tslen=num_timepoints, highpass_freq=highpass_freq, target_ta=max(0,ta_delta1)) for ta_delta1 in ta_delta1s]
67. # Use these alpha values to construct desired frequency spectra
68. spectra = np.asarray([make_spectrum(tslen=num_timepoints, sample_rate=sample_rate, alpha=alpha, highpass_freq=highpass_freq) for alpha in alphas])
69. # Spatial autocorrelation matrix
70. corr = spatial_exponential_floor(distance_matrix, sa_lambda, sa_inf)
71. # Compute timeseries from desired correlation matrix and frequency spectra
72. tss = spatial_temporal_timeseries(cm=corr, spectra=spectra, seed=seed)
73. return tss
76. def spatial_temporal_timeseries(cm, spectra, seed=None):
77. """Generate timeseries with given amplitude spectra and correlation matrices
79. Args:
80. cm (NxN numpy array): The correlation matrix
81. spectra (Nxk numpy array): A list of Fourier spectra generated by [make_spectrum][spatiotemporal.models.make_spectrum]
82. Each of the N spectra are associated with a row/column of `cm`.
83. seed (int, optional): the random seed. If not specified, it will use
84. the current state of the numpy random number generator.
86. Returns:
87. NxT numpy array: N timeseries of length T. Timeseries i will have a power spectrum given by `spectra[i]`,
88. and will be correlated with the other timeseries with correlations cm[i].
89. """
90. N_regions = cm.shape[0]
91. N_freqs = len(spectra[0])
92. N_timepoints = (N_freqs-1)*2
93. assert spectra.shape == (N_regions, N_freqs)
94. sum_squares = np.sum(spectra**2, axis=1, keepdims=True)
95. cosine_similarity = (spectra @ spectra.T)/np.sqrt(sum_squares @ sum_squares.T)
96. covmat = cm / cosine_similarity
97. if np.min(np.linalg.eigvalsh(covmat)) < -1e-8:
98. raise PositiveSemidefiniteError("Correlation matrix is not possible with those spectra using this method!")
99. randstate = np.random.RandomState(seed)
100. rvs = randstate.multivariate_normal(np.zeros(N_regions), cov=covmat, size=N_freqs*2)
101. reals = rvs[0:N_freqs].T * spectra
102. ims = rvs[N_freqs:].T * spectra
103. # Since the signal length is even, frequencies +/- 0.5 are equal
104. # so the coefficient must be real.
105. ims[:,-1] = 0
106. # The DC component must be zero and real
107. reals[:,0] = 0
108. ims[:,0] = 0
109. tss = np.fft.irfft(reals + 1J*ims, n=N_timepoints, axis=-1)
110. return tss
112. def make_spectrum(tslen, sample_rate, alpha, highpass_freq):
113. """Create a 1/f^alpha spectrum.
115. Args:
116. tslen (int): the length of the timeseries represented by the spectrum
117. sample_rate (float): the spacing between timepoints (e.g. TR in fMRI)
118. alpha (float): the pink noise exponent, between 0 and 2.
119. highpass_freq (float): if non-zero, apply a highpass filter above the
120. given frequency. A good default is 0.01 for fMRI timeseries.
122. Returns:
123. length-k numpy array: Return the fourier spectrum (amplitude spectrum)
124. """
125. freqs = np.fft.rfftfreq(tslen, sample_rate)
126. with np.errstate(all="warn"):
127. spectrum = freqs**(-alpha/2)
128. if highpass_freq > 0:
129. butter = scipy.signal.iirfilter(ftype='butter', N=4, Wn=highpass_freq, btype='highpass', fs=1/sample_rate, output='ba')
130. butterresp = scipy.signal.freqz(*butter, fs=1/sample_rate, worN=len(freqs), include_nyquist=True)
131. assert np.all(np.isclose(freqs, butterresp[0]))
132. spectrum = spectrum * np.abs(butterresp[1])
133. spectrum[0] = 0
134. return spectrum
136. def how_much_noise(spectrum, target_ta):
137. """Determine the standard deviation of noise to add to achieve a target TA-Δ₁.
139. This function answers the following question: If I generate timeseries with
140. frequency spectrum (amplitude spectrum) `spectrum`, and then add white
141. noise to the generated timeseries, what should the standard deviation of
142. this white noise be if I want the timeseries to have the TA-Δ₁ coefficient
143. `target_ta`?
145. Args:
146. spectrum (length-k numpy array): the power spectrum to generate from, e.g.,
147. that generated from [make_spectrum][spatiotemporal.models.make_spectrum].
148. target_ta (float): the desired
149. TA-Δ₁.
151. Returns:
152. float: The standard deviation of white noise to add
153. """
155. N = len(spectrum)
156. weightedsum = np.sum(spectrum[1:]**2*np.cos(np.pi*np.arange(1, N)/N))
157. try:
158. sigma = np.sqrt((weightedsum - np.sum(spectrum**2)*target_ta)/(target_ta*N**2))
159. except FloatingPointError:
160. sigma = 0
161. return sigma
163. def ta_to_alpha(tslen, sample_rate, highpass_freq, target_ta):
164. """Compute the (pink noise) alpha which would give, noiseless, the given TA-Δ₁.
166. Generate timeseries with get_spectrum_ta, i.e. high pass filtered.
168. Args:
169. tslen (int): the length of the timeseries represented by the spectrum
170. sample_rate (float): the spacing between timepoints (e.g. TR in fMRI)
171. alpha (float): the pink noise exponent, between 0 and 2.
172. highpass_freq (float): if non-zero, apply a highpass filter above the
173. given frequency. A good default is 0.01 for fMRI timeseries.
175. Returns:
176. float: a value of alpha such that the filtered pink noise with
177. this exponent has TA-Δ₁ coefficient `target_ta`.
178. """
179. objfunc = lambda alpha : (get_spectrum_ta(make_spectrum(tslen, sample_rate, alpha[0], highpass_freq)) - target_ta)**2
180. x = scipy.optimize.minimize(objfunc, 1.5, bounds=[(0, 2)])
181. return float(x.x[0])
183. def get_spectrum_ta(spectrum):
184. """Given a fourier spectrum, return the expected TA-Δ₁ of a timeseries generated with that spectrum.
186. Args:
187. spectrum (length-k numpy array): the power spectrum to generate from, e.g.,
188. that generated from [make_spectrum][spatiotemporal.models.make_spectrum].
190. Returns:
191. float: the TA-Δ₁ value that would be expected if a timeseries had the given power spectrum and random phases.
193. """
194. N = len(spectrum)
195. weightedsum = np.sum(spectrum**2*np.cos(np.pi*np.arange(0, N)/N))
196. return weightedsum/np.sum(spectrum**2)
198. ta_to_alpha_cache = {}
199. def ta_to_alpha_fast(tslen, sample_rate, highpass_freq, target_ta):
200. """Identical to `ta_to_alpha`, but discretize and cache to increase speed.
202. See [ta_to_alpha][spatiotemporal.models.ta_to_alpha] for documentation.
203. """
204. global ta_to_alpha_cache
205. taround = round(target_ta, 2)
206. key = (tslen, sample_rate, highpass_freq, taround)
207. if key in ta_to_alpha_cache.keys():
208. return ta_to_alpha_cache[key]
209. val = ta_to_alpha(*key)
210. ta_to_alpha_cache[key] = val
211. return val
213. class PositiveSemidefiniteError(Exception):
214. pass

models.py at commit c00dec6, no license · at the source