**Reference**

HowToGuides

Manuals

LabAlumni

DataAnalysis

Advice for...

**Admin**

analysis:course-w16:week10

Goals:

- Learn to estimate and plot tuning curves, raw and smoothed
- Implement a basic Bayesian decoding algorithm
- Compare decoded and actual position by computing the decoding error

Resources:

- (if you have not encountered Bayes' rule and conditional probability before) A brief introduction by Ken Murphy at UBC, and one from the Scientific American
- (for reference) Zhang et al. 1998, first application of decoding to place cell data, with nice explanations and derivations
- (for reference) Brown et al. 1998, an example of a more sophisticated decoding method

To support adaptive behavior, activity in the brain must correspond in some way to relevant sensory events and planned movements, combine many sources of information into multimodal percepts, and recall traces of past events to inform predictions about the future. In other words, neural activity must somehow *encode* relevant quantities. For instance, it can be demonstrated behaviorally that many animals use estimates of their location and head direction to navigate towards a goal. Where, and how, are these quantities represented in the brain? What are the neural circuits that can compute and update these signals? How do place and direction estimates contribute to which way to go?

This *information processing* view of the brain has been extremely influential, as highlighted by the enduring appeal of Hubel and Wiesel's demonstrations that single cells in macaque V1 respond to bars of light not only within a particular region of visual space, but also with a specific orientation. Such cells are said to be *tuned* for orientation [of the bar] and a typical *tuning curve* would therefore look like this:

This tuning curve describes how the cell responds, on average, to different orientations of the stimulus. If the cell were to respond with the same firing rate across the range of stimulus orientations, then the cell is indifferent to this particular stimulus dimension: it does not encode it. However, because this cell clearly modulates its firing rate with stimulus orientation, it encodes, or represents (I use these terms interchangeably, but some disagree) this quantity in its activity.

We can turn this idea around and note that if orientation is encoded, this implies we can also *decode* the original stimulus from the cell's activity. For instance, if we noted that this cell was firing at a high rate, we would infer that the stimulus orientation is likely close to the cell's preferred direction. Note that this requires knowledge of the cell's tuning curve, and that based on one cell only, we are unlikely to be able to decode (or reconstruct, which means the same thing) the stimulus perfectly. The more general view is to say that the cell's activity provides a certain amount of information about the stimulus, or equivalently, that our (decoded) estimate of the stimulus is improved by taking the activity of this cell into account.

This module first explores some practical issues in estimating tuning curves of “place cells” recorded from the rat hippocampus. An introduction to a particular decoding method (Bayesian decoding) is followed by application to many simultaneously recorded place cells as a rat performs a T-maze task.

First, load the “place cell” data set also used in the previous module, which contains a number of spike trains recorded simultaneously from the dorsal CA1 area of the hippocampus:

%% cd('D:\data\R042\R042-2013-08-18'); please = []; please.load_questionable_cells = 1; S = LoadSpikes(please); pos = LoadPos([]);

The `load_questionable_cells`

option in `LoadSpikes()`

results in the loading of `*._t`

files, in addition to the familiar `*.t`

spike time files. The underscore extension indicates a cell with questionable isolation quality, likely contaminated with noise, spikes from other neurons, and/or missing spikes. In general, you do not want to use such neurons for analysis, but in this case we are not concerned with properties of individual neurons. We are instead interested in the information present in a population of neurons, and for this we will take everything we can get!

Before looking at the data, we will first exclude the pre- and post-run segments of the data:

LoadExpKeys; S = restrict(S,ExpKeys.TimeOnTrack,ExpKeys.TimeOffTrack); pos = restrict(pos,ExpKeys.TimeOnTrack,ExpKeys.TimeOffTrack);

Now we can plot the position data:

plot(getd(pos,'x'),getd(pos,'y'),'.','Color',[0.5 0.5 0.5],'MarkerSize',1); axis off; hold on;

Note that `getd()`

is a utility function that retrieves data associated with a specific label; see Module 2 for details.

Next, we plot the spikes of a single cell at the location where the rat was when each spike was emitted:

iC = 7; spk_x = interp1(pos.tvec,getd(pos,'x'),S.t{iC},'linear'); spk_y = interp1(pos.tvec,getd(pos,'y'),S.t{iC},'linear'); h = plot(spk_x,spk_y,'.r');

Note the use of `interp1()`

here: it finds the corresponding `x`

and `y`

values for each spike time using linear interpolation. You should see:

This cell seems to have a place field just to the left of the choice point on the track (a T-maze). There are also a few spikes on the pedestals, where the rat rests in between runs on the track.

This figure is a useful visualization of the raw data, but it is not a tuning curve. As a first step towards estimating this cell's tuning curve (or *encoding model*, we should restrict the spikes to only those occurring when the rat is running on the track:

LoadMetadata; ENC_S = restrict(S,metadata.taskvars.trial_iv); ENC_pos = restrict(pos,metadata.taskvars.trial_iv); % check for empties and remove keep = ~cellfun(@isempty,ENC_S.t); ENC_S.t = ENC_S.t(keep); ENC_S.label = ENC_S.label(keep); S.t = S.t(keep); S.label = S.label(keep);

We have created `ENC_`

versions of our spike trains and position data, containing only data from when the rat was running on the track (using experimenter annotation stored in the metadata; `trial_iv`

contains the start and end times of trials) and removed all cells from the data set that did not have any spikes on the track.

☛ Plot the above scatterfield again for the restricted spike train. Verify that no spikes are occurring off the track by comparing your plot to the previous one for the full spike trains, above.

To estimate tuning curves from the data, we need to divide spike count by time spent for each location on the maze. A simple way of doing that is to obtain 2-D histograms, shown here for the position data:

clear pos_mat; pos_mat(:,1) = getd(ENC_pos,'y'); % construct input to 2-d histogram pos_mat(:,2) = getd(ENC_pos,'x'); SET_xmin = 80; SET_ymin = 0; % set up bins SET_xmax = 660; SET_ymax = 520; SET_xBinSz = 10; SET_yBinSz = 10; x_edges = SET_xmin:SET_xBinSz:SET_xmax; y_edges = SET_ymin:SET_yBinSz:SET_ymax; occ_hist = histcn(pos_mat,y_edges,x_edges); % 2-D version of histc() no_occ_idx = find(occ_hist == 0); % NaN out bins never visited occ_hist(no_occ_idx) = NaN; occ_hist = occ_hist .* (1/30); % convert samples to seconds using video frame rate (30 Hz) subplot(221); pcolor(occ_hist); shading flat; axis off; colorbar title('occupancy');

We can do the same thing for the spikes of our example neuron:

% basic spike histogram clear spk_mat; iC = 7; spk_x = interp1(ENC_pos.tvec,getd(ENC_pos,'x'),ENC_S.t{iC},'linear'); spk_y = interp1(ENC_pos.tvec,getd(ENC_pos,'y'),ENC_S.t{iC},'linear'); spk_mat(:,2) = spk_x; spk_mat(:,1) = spk_y; spk_hist = histcn(spk_mat,y_edges,x_edges); spk_hist(no_occ_idx) = NaN; subplot(222) pcolor(spk_hist); shading flat; axis off; colorbar title('spikes');

..and finally simply divide one by the other:

% rate map tc = spk_hist./occ_hist; subplot(223) pcolor(tc); shading flat; axis off; colorbar title('rate map');

This gives:

Note that from the occupancy map, you can see the rat spent relatively more time at the base of the stem compared to other segments of the track. However, the rough binning is not very satisfying. Let's see if we can do better with some smoothing:

kernel = gausskernel([4 4],2); % Gaussian kernel of 4x4 pixels, SD of 2 pixels (note this should sum to 1) [occ_hist,~,~,pos_idx] = histcn(pos_mat,y_edges,x_edges); occ_hist = conv2(occ_hist,kernel,'same'); occ_hist(no_occ_idx) = NaN; occ_hist = occ_hist .* (1/30); subplot(221); pcolor(occ_hist); shading flat; axis off; colorbar title('occupancy'); % spk_hist = histcn(spk_mat,y_edges,x_edges); spk_hist = conv2(spk_hist,kernel,'same'); % 2-D convolution spk_hist(no_occ_idx) = NaN; subplot(222) pcolor(spk_hist); shading flat; axis off; colorbar title('spikes'); % tc = spk_hist./occ_hist; subplot(223) pcolor(tc); shading flat; axis off; colorbar title('rate map');

Now you should get:

These are well-formed tuning curves we can use for decoding. Of course we could bin more finely for increased spatial resolution, but this will slow down the decoding, so for now it's not worth it.

Next we obtain a tuning curve for all our cells:

clear tc all_tc nCells = length(ENC_S.t); for iC = 1:nCells spk_x = interp1(ENC_pos.tvec,getd(ENC_pos,'x'),ENC_S.t{iC},'linear'); spk_y = interp1(ENC_pos.tvec,getd(ENC_pos,'y'),ENC_S.t{iC},'linear'); clear spk_mat; spk_mat(:,2) = spk_x; spk_mat(:,1) = spk_y; spk_hist = histcn(spk_mat,y_edges,x_edges); spk_hist = conv2(spk_hist,kernel,'same'); spk_hist(no_occ_idx) = NaN; tc = spk_hist./occ_hist; all_tc{iC} = tc; end

We can inspect the results as follows:

%% ppf = 25; % plots per figure for iC = 1:length(ENC_S.t) nFigure = ceil(iC/ppf); figure(nFigure); subtightplot(5,5,iC-(nFigure-1)*ppf); pcolor(all_tc{iC}); shading flat; axis off; caxis([0 10]); end

You will see a some textbook “place cells” with a clearly defined single place field. There are also cells with other firing patterns.

The data in this module computes tuning curves for location, but the idea is of course more general. For continuous variables in particular, it is a natural and powerful way to describe the relationship between two quantities – spikes and location in this case, but there is no reason why you couldn't do something like pupil diameter as a function of arm reaching direction, for instance!

As noted in the introduction above, given that we have neurons whose activity seems to encode some stimulus variable (location in this case), we can attempt to decode that variable based on the neurons' time-varying activity.

A popular approach to doing this is “one-step Bayesian decoding”, illustrated in this figure (from van der Meer et al. 2010):

For this particular experiment, the goal of decoding is to recover the location of the rat, given neural activity in some time window. More formally, we wish to know $P(\mathbf{x}|\mathbf{n})$, the probability of the rat being at each possible location $x_i$ ($\mathbf{x}$ in vector notation, to indicate that there are many possible locations) given a vector of spike counts $\mathbf{n}$.

If $P(\mathbf{x}|\mathbf{n})$ (the “posterior”) is the same for every location bin $x_i$ (i.e. is uniform), that means all locations are equally likely and we don't have a good guess; in contrast, if most of the $x_i$ are zero and a small number have a high probability, that means we are confident predicting the most likely location. Of course, there is no guarantee that our decoded estimate will agree with the actual location; we will test this later on.

So how can we obtain $P(\mathbf{x}|\mathbf{n})$? We can start with Bayes' rule:

\[P(\mathbf{x}|\mathbf{n})P(\mathbf{n}) = P(\mathbf{n}|\mathbf{x})P(\mathbf{x})\]

If you have not come across Bayes' rule before, or the above equation looks mysterious to you, review the gentle intro by linked to at the top of the page. In general, it provides a quantitative way to update prior beliefs in the face of new evidence.

The key quantity to estimate is $P(\mathbf{n}|\mathbf{x})$, the probability of observing $n$ spikes in a given time window when the rat is at location $x$. At the basis of estimating this probability (the “likelihood” or evidence) lies the tuning curve: this tells us the *average* firing rate at each location. We need a way to convert a given number of spikes – whatever we observe in the current time window for which we are trying to decode activity, 3 spikes for cell 1 in the figure above – to a probability. In other words, what is the probability of observing 3 spikes in a 250ms time window, given that for this location the cell fires, say at 5Hz on average?

A convenient answer is to assume that the spike counts follow a Poisson distribution. Assuming this enables us to assign a probability to each possible spike count for a mean firing rate given by the tuning curve. For instance, here are the probabilities of observing different numbers of spikes $k$ (on the horizontal axis) for four different means ($\lambda = $1, 4 and 10):

In general, from the definition of the Poisson distribution, it follows that

\[P(n_i|\mathbf{x}) = \frac{(\tau f_i(\mathbf{x}))^{n_i}}{n_i!} e^{-\tau f_i (\mathbf{x})}\]

$f_i(\mathbf{x})$ is the average firing rate of neuron $i$ over $x$ (i.e. the tuning curve for position), $n_i$ is the number of spikes emitted by neuron $i$ in the current time window, and $\tau$ is the size of the time window used. Thus, $\tau f_i(\mathbf{x})$ is the mean number of spikes we expect from neuron $i$ in a window of size $\tau$; the Poisson distribution describes how likely it is that we observe the actual number of spikes $n_i$ given this expectation.

In reality, place cell spike counts are typically not Poisson-distributed ( Fenton et al. 1998) so this is clearly a simplifying assumption. There are many other, more sophisticated approaches for the estimation of $P(n_i|\mathbf{x})$ (see for instance Paninski et al. 2007) but this basic method works well for many applications.

The above equation gives the probability of observing $n$ spikes for a given average firing rate for a single neuron. How can we combine information across neurons? Again we take the simplest possible approach and assume that the spike count probabilities for different neurons are independent. This allows us to simply multiply the probabilities together to give:

\[P(\mathbf{n}|\mathbf{x}) = \prod_{i = 1}^{N} \frac{(\tau f_i(\mathbf{x}))^{n_i}}{n_i!} e^{-\tau f_i (\mathbf{x})}\]

An analogy here is simply to ask: if the probability of a coin coming up heads is $0.5$, what is the probability of two coints, flipped simultaneously, coming up heads? If the coins are independent then this is simply $0.5*0.5$.

Combining the above with Bayes' rule, and rearranging a bit, gives

\[P(\mathbf{x}|\mathbf{n}) = C(\tau,\mathbf{n}) P(\mathbf{x}) (\prod_{i = 1}^{N} f_i(\mathbf{x})^{n_i}) \: e (-\tau \sum_{i = 1}^N f_i(\mathbf{x})) \]

This is more easily evaluated in vectorized MATLAB code. $C(\tau,\mathbf{n})$ is a normalization factor which we simply set to guarantee $\sum_x P(\mathbf{x}|\mathbf{n}) = 1$ (Zhang et al. 1998). For now, we assume that $P(\mathbf{x})$ (the “prior”) is uniform, that is, we have no prior information about the location of the rat and let our estimate be completely determined by the likelihood.

The tuning curves take care of the $f_i(x)$ term in the decoding equations. Next, we need to get $\mathbf{n}$, the spike counts.

To obtain spike counts within a given bin size, we can use histc():

%% make Q-mat clear Q; binsize = 0.25; % seconds % assemble tvecs tvec_edges = metadata.taskvars.trial_iv.tstart(1):binsize:metadata.taskvars.trial_iv.tend(end); Q_tvec_centers = tvec_edges(1:end-1)+binsize/2; for iC = length(ENC_S.t):-1:1 spk_t = ENC_S.t{iC}; Q(iC,:) = histc(spk_t,tvec_edges); Q(iC,end-1) = Q(iC,end-1)+Q(iC,end); % remember last bin of histc() gotcha end Q = Q(:,1:end-1);

This “Q-matrix” of size `[nCells x nTimeBins]`

is the start of a number of analyses, such as the nice ensemble reactivation procedure introduced in Peyrache et al. 2009. Let's inspect it briefly:

imagesc(Q_tvec_centers,1:nCells,Q) set(gca,'FontSize',16); xlabel('time(s)'); ylabel('cell #');

Our Q-matrix only includes non-zero counts when the animal is running on the track; these episodes manifest as narrow vertical stripes. To speed up calculations later, let's restrict Q to those times only:

%% only include data we care about (runs on the maze) Q_tsd = tsd(Q_tvec_centers,Q); Q_tsd = restrict(Q_tsd,metadata.taskvars.trial_iv);

The final step before the actual decoding procedure is to reformat the tuning curves a bit to make the decoding easier to run. Instead of keeping them as a 2-D matrix, we just unwrap this into 1-D:

%% prepare tuning curves clear tc nBins = numel(occ_hist); nCells = length(S.t); for iC = nCells:-1:1 tc(:,:,iC) = all_tc{iC}; end tc = reshape(tc,[size(tc,1)*size(tc,2) size(tc,3)]); occUniform = repmat(1/nBins,[nBins 1]);

Aaandd… action!

%% decode Q_tvec_centers = Q_tsd.tvec; Q = Q_tsd.data; nActiveNeurons = sum(Q > 0); len = length(Q_tvec_centers); p = nan(length(Q_tvec_centers),nBins); for iB = 1:nBins tempProd = nansum(log(repmat(tc(iB,:)',1,len).^Q)); tempSum = exp(-binsize*nansum(tc(iB,:),2)); p(:,iB) = exp(tempProd)*tempSum*occUniform(iB); end p = p./repmat(sum(p,2),1,nBins); % renormalize to 1 total probability p(nActiveNeurons < 1,:) = 0; % ignore bins with no activity

☛ Compare these steps with the equations above.

- There is no log in the equations; why does it appear here?
- What parts of the equation correspond to the
`tempProd`

and`tempSum`

variables?

The hard work is done. Now we just need to display the results. Before we do so, we should convert the rat's actual position into our binned form, so that we can compare it to the decoded estimate:

xBinned = interp1(ENC_pos.tvec,pos_idx(:,1),Q_tvec_centers); yBinned = interp1(ENC_pos.tvec,pos_idx(:,2),Q_tvec_centers);

Now we can visualize the decoding (press Ctrl-C to break out of the loop):

goodOccInd = find(occ_hist > 0); SET_nxBins = length(x_edges)-1; SET_nyBins = length(y_edges)-1; dec_err = nan(length(Q_tvec_centers),1); for iT = 1:length(Q_tvec_centers) cla; temp = reshape(p(iT,:),[SET_nyBins SET_nxBins]); toPlot = nan(SET_nyBins,SET_nxBins); toPlot(goodOccInd) = temp(goodOccInd); pcolor(toPlot); axis xy; hold on; caxis([0 0.5]); shading flat; axis off; hold on; plot(yBinned(iT),xBinned(iT),'ow','MarkerSize',15); % get x and y coordinates of MAP [~,idx] = max(toPlot(:)); [x_map,y_map] = ind2sub(size(toPlot),idx); if nActiveNeurons(iT) > 0 dec_err(iT) = sqrt((yBinned(iT)-y_map).^2+(xBinned(iT)-x_map).^2); end plot(y_map,x_map,'g*','MarkerSize',5); h = title(sprintf('t %.2f, nCells %d, dist %.2f',Q_tvec_centers(iT),nActiveNeurons(iT),dec_err(iT))); if nActiveNeurons(iT) == 0 set(h,'Color',[1 0 0]); else set(h,'Color',[0 0 0]); end drawnow; pause(0.1); end

This plot shows the posterior $P(\mathbf{x}|\mathbf{n})$, as the rat moves around the maze; its actual position is indicated by the white `o`

, and the pixel with the highest posterior probability is indicated by the green `*`

. As you can see, the decoding seems to track the rat's actual location as it moves.

☛ No decoding is available for those bins where no neurons are active, because we manually set the posterior to zero. However, there also seem to be some frames in the animation where some neurons are active (as indicated in the title), yet no decoded estimate is visible. What is the explanation for this?

By making the MATLAB animation into a movie file, it is often easier to explore the results. To do this, we can run the animation code above, with a few small modifications. First, before entering the main plotting loop, set the figure to be used to a specific size:

h = figure; set(h,'Position',[100 100 320 240]);

This is important first, to keep the size of the resulting movie file manageable (the above sets a 320×240 pixel figure size), and second, because many movie encoders (such as the excellent XVid) will only work with certain sizes.

Next, we need to store each frame into a variable that we can later write to file. Modify the last two lines inside the loop to:

f(iT) = getframe(gcf); % store current frame drawnow;

If you now run the code again, each frame gets stored in the `f`

variable as the loop runs. Break out of the loop after a few seconds to test the writing-to-file part:

fname = 'test.avi'; movie2avi(f,fname,'COMPRESSION','XVid','FPS',10);

The above will only work if you have the XVid codec installed: I highly recommend this because it creates movie files that are an order of magnitude smaller than uncompressed files. If you have trouble with XVid, you can of course still save an uncompressed file for now. For longer movies, it is often required to save a file, say, every 500 frames, to prevent the `f`

variable getting too large. These segments can then be merged with a video editing program such as VirtualDub (Windows only AFAIK; please suggest OSX/Linux alternatives if you know any that work well!).

By running the above loop without plotting we can obtain the “decoding error”, that is the distance between the maximum a posteriori (MAP) estimate and the rat's true position:

%% get distance (no plotting) dec_err = nan(length(Q_tvec_centers),1); SET_nxBins = length(x_edges)-1; SET_nyBins = length(y_edges)-1; xBinned = interp1(ENC_pos.tvec,pos_idx(:,1),Q_tvec_centers); yBinned = interp1(ENC_pos.tvec,pos_idx(:,2),Q_tvec_centers); for iT = 1:length(Q_tvec_centers); temp = reshape(p(iT,:),[SET_nyBins SET_nxBins]); toPlot = nan(SET_nyBins,SET_nxBins); toPlot(goodOccInd) = temp(goodOccInd); % get x and y coordinates of MAP [~,idx] = max(toPlot(:)); [x_map,y_map] = ind2sub(size(toPlot),idx); if nActiveNeurons(iT) > 0 dec_err(iT) = sqrt((yBinned(iT)-y_map).^2+(xBinned(iT)-x_map).^2); end end

A nice way to plot this is to average by lap as well as overall:

% get trial id for each sample trial_id = zeros(size(Q_tvec_centers)); trial_idx = nearest_idx3(metadata.taskvars.trial_iv.tstart,Q_tvec_centers); % NOTE: on non-Windows, use nearest_idx.m trial_id(trial_idx) = 1; trial_id = cumsum(trial_id); %plot(dec_err_tsd.tvec,trial_id,'.k'); figure; set(gca,'FontSize',18); boxplot(dec_err,trial_id); xlabel('trial'); ylabel('decoding error (pixels)'); av_error = nanmean(dec_err); title(sprintf('avg err %.2f',av_error));

This yields:

(Note, your plot might look a little different.)

Thus, on average our estimate is 2.14 pixels away from the true position. Earlier laps seem to have some more outliers of bins where our estimate is bad (large distance) but there is no obvious trend across laps visible.

☛ How does the decoding accuracy depend on the bin size used? Try a range from very small (10ms) to very large (1s) bins, making sure to note the average decoding error for 50ms bins, for comparison with results in the next module. What factors need to be balanced if the goal is maximum accuracy?

A different way of looking at the decoding error is to plot it as a function of space:

cfg = []; cfg.y_edges = y_edges; cfg.x_edges = x_edges; dec_err_tsd = tsd(Q_tvec_centers,dec_err); space_err = TSDbySpace(cfg,ENC_pos,dec_err_tsd); figure; pcolor(space_err); shading flat; axis off; colorbar; caxis([0 10]);

This gives:

It looks like the decoding error is on average larger on the central stem, compared to the arms of the maze.

☛ What could be some reasons for this? Can you think of ways to test your suggestions?

★ Implement a decoding analysis on your own data. Remember that this does not necessarily requires using spiking data – anything that you can construct a tuning curve for would work! In this module, we had something like 100 simultaneously recorded neurons, but even if you have only one, you can still attempt to use it for decoding. Quantify decoding performance (error) for a few relevant parameters.

★ How does decoding performance scale with the number of cells used? This is an important issue if we want to figure out if we should invest resources in attempting to record from more neurons, or if we have all we need in data sets such as this one.

★ In a famous paper, Johnson and Redish (2007) showed that the hippocampus transiently represents possible future trajectories as rats appeared to deliberate between choices (left? right?) at a decision point. However, they used a controversial “two-step” decoding algorithm which attracted criticism. Refer to the Methods section of that paper to figure out how they did the decoding, and modify the code above to implement their version. What differences do you notice?

analysis/course-w16/week10.txt · Last modified: 2018/07/07 10:19 (external edit)

Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 4.0 International

## Discussion

What is the variable run_start in the quantifying decoding accuracy section?

oops, that shouldn't have been there – fixed! (it was a temporary variable containing trial start times..)

generic tadalafil https://nextadalafil.com/ tadalis sx <a href=“https://nextadalafil.com/#”>tadalafil daily use</a>

generic tadalafil 40 mg https://pulmoprestadalafil.com/ buy tadalafil us <a href=“https://pulmoprestadalafil.com/#”>tadalafil tablets</a>

buy tadalafil https://extratadalafill.com/ tadalafil 30 mg <a href=“https://extratadalafill.com/#”>tadalafil generic</a>

tadalafil 30 mg https://superactivetadalafil.com/ tadalafil 60 mg for sale <a href=“https://superactivetadalafil.com/#”>tadalafil 30 mg</a>

40 mg tadalafil https://tadalafilgenc.com/ tadalis sx <a href=“https://tadalafilgenc.com/#”>tadalis sx</a>

buy tadalafil https://tadalafilgenc.com/ tadalafil gel <a href=“https://tadalafilgenc.com/#”>tadalafil 40</a>

tadalafil max dose https://pulmoprestadalafil.com/ tadalafil gel <a href=“https://pulmoprestadalafil.com/#”>tadalafil max dose</a>

generic tadalafil united states https://extratadalafill.com/ tadalafil 40 <a href=“https://extratadalafill.com/#”>tadalafil pills 20mg</a>

tadalafil generic https://superactivetadalafil.com/ generic tadalafil united states <a href=“https://superactivetadalafil.com/#”>40 mg tadalafil</a>

generic tadalafil 40 mg https://nextadalafil.com/ tadalafil max dose <a href=“https://nextadalafil.com/#”>tadalafil 30 mg</a>

href=“https://tanphuchung.vn/dieu-hoa-multi-la-gi-don-vi-lap-dat-dieu-hoa-multi-uy-tin/”>điều hòa multi</a>, <a href=“https://tanphuchung.vn/san-pham/dieu-hoa-daikin-vrv/”>điều hòa vrv</a>, <a href=“https://tanphuchung.vn/dich-vu-thiet-ke-va-thi-cong-he-thong-dieu-hoa-trung-tam-vrv/”>hệ thống điều hòa trung tâm vrv</a>, <a href=“https://tanphuchung.vn/quy-trinh-tu-van-lap-dat-dieu-hoa-trung-tam-vrv-cho-biet-thu/”>lắp đặt điều hòa vrv</a> <a href=“https://tanphuchung.vn/dich-vu-lap-dat-dieu-hoa-vrv-nhanh-tai-ha-noi-tan-phuc-hung/”>lắp đặt điều hòa vrv tại hà nội</a> <a href=“https://tanphuchung.vn/cach-lap-dieu-hoa-multi-cho-chung-cu-3-phong-ngu-1-phong-khach/”>điều hòa multi cho chung cư</a> <a href=“https://tanphuchung.vn/quy-trinh-tu-van-lap-dat-dieu-hoa-trung-tam-vrv-cho-biet-thu/”>điều hòa trung tâm biệt thự</a> <a href=“https://tanphuchung.vn/quy-trinh-tu-van-thiet-ke-thi-cong-dieu-hoa-multi-cho-biet-thu/”>điều hòa multi cho biệt thự</a> <a href=“https://tanphuchung.vn/san-pham/dieu-hoa-dan-dung-multi-mitsubishi-heavy/”>điều hòa multi mitsubishi</a> <a href=“https://tanphuchung.vn/san-pham/dieu-hoa-trung-tam-vrf-mitsubishi-heavy/”>điều hòa vrv mitsubishi</a>

KHO LẠNH ĐẠI AN CHUYÊN THI CÔNG LẮP ĐẶT KHO LẠNH Địa chỉ: Số 16, ngõ 103, đường Phương Canh, Phường Xuân Phương, Nam Từ Liêm, Hà Nội Điện Thoại : 0982.634.000 Hotline: 0982.634.000 Website: https://kholanhdaian.com Email:kholanhdaian@gmail.com <p><b>xem thêm: => </b><a href=“https://kholanhdaian.com/bang-bao-gia-lap-dat-kho-lanh-tot-nhat-thi-truong-2021/”><b>bảng giá kho lạnh</b></a></p> <p><b>xem thêm: => </b><a href=“https://kholanhdaian.com/kho-lanh-bao-quan-la-gi-don-vi-lap-dat-kho-lanh-uy-tin/”><b>kho lạnh</b></a></p> <p><b>xem thêm: => </b><a href=“https://kholanhdaian.com/kho-lanh-bao-quan-la-gi-don-vi-lap-dat-kho-lanh-uy-tin/”><b>lắp đặt kho lạnh</b></a></p>