Version from 6.10.2014

Command reference: Matlab syntax

pdf version of the script: [Tag 8.pdf]

Word version of the script: [Tag 8.docx]

Downloads:

T8A6) [stimT8.mat] [ratenT8.mat]

T8A7) [allereizeT8.mat] [allrateT8.mat]

T8H1) [populationen.mat]

T8H2) [intfire.m]

Topics:

Curve fitting is used to describe measurement data using a curve - i.e. a mathematically describable function. The idea is that the measurement data is systematically dependent on a given parameter, but is statistically adjusted by the expected values.

• The values of the parameter (e.g. time, fertiliser quantity, sunshine duration, etc.) are plotted graphically on the x-axis, the measured values and the theoretical values as a function of this parameter on the y-axis.
• The measured values are often shown as dots (or other symbols) and the theoretical values as a solid curve.
• This representation helps to understand the idea of interpolation: It can be used to estimate what measured values would have been expected if additional measurement points had been determined between the selected parameter values.
• In addition, curve fitting also enables extrapolation, the prediction of how the measured values will behave outside the tested parameter range.

In order to find a function that describes the dependence of the existing measurement data on the varied parameter as well as possible, you want to find a curve that is as close as possible to the measurement data. For this (or generally when trying to describe measurement data using a formation rule (i.e. a model)), you need a measure of how well the model reflects the data. This allows you to compare different models with each other and choose the one that fits best. A standard procedure to determine the optimal function (or model) is to minimise the mean square deviation between the measured values and the estimated values. (In principle, however, there are many other ways to measure the quality of a model.)

To do this, the residual is first determined for each data point, i.e. the difference between the empirically determined value y i and the estimate for this value (estimated values are often labelled with a roof):

The mean squared deviation between all data points and the estimate is determined as the sum of the squared residuals, divided by the number of data points n:

• What does the square do?

In this course we use polynomials of the form

to fit a function to the empirical data.

• What does a straight line with slope 2.5 and y-intercept 7 look like in this illustration?
• Of course, not all dependencies can be adequately represented by polynomials. What counterexamples can you think of?

Matlab: In Matlab, the commands polyfit and polyval are available for curve fitting of a polynomial to measurement data: [p,S]=polyfit(x,y,n)

• Return value p=[p1 p2... pn pn+1] returns the coefficients for the equation

y_schaetz = polyval(p,x,S) receives as input

• the value given by polyfit returns the vector p
• returned by polyfit returned structure S,
• and a vector of x-values.

For thesexvaluespolyvalcalculates the y values predicted according to the equation found (the x values are inserted into the equation). [ y_schaetz,delta] = polyval(p,x,S) also provides the output argument

• delta as a measure of how reliable this estimate is. In the range y_schaetz-delta to y_schaetz+delta are 50% of the measuring points.

Tasks: I would like to ask you to first carry out curve fitting "by hand" for a simple case so that you can better understand the basic idea: T8A1) In preparation, let's first deal with the error measure, the mean square error: write a function that gets a vector of measured y-values and a vector of estimated y-values of the same length and returns the mean square error. *T8A2) Extend your function so that it also calculates the mean square error for a matrix of measured values (where each x-value is assigned several y-values) and a vector of estimated y-values (one for each x-value). T8A3) Now for the actual curve fitting:
When nerve cells are stimulated with electricity, they respond with a sequence of action potentials.
Imagine a neuron that reacts to a current pulse of 1 nA with a frequency of 32 sp/s, at 2 nA it produces 60 sp/s.

• What response frequency do you expect with a current injection of 1.5 nA?
• What is the linear equation for the straight line connecting the two points?
• Plot the straight line and the two data points

T8A4) You now carry out another measurement and realise that 44 sp/s are generated at 1.5 nA current injection.

• Expand your plot to include the new data point.
• Calculate the square vertical distance of the data point from the straight line.
• Estimate "by eye" a straight line that approximates all three points well and plot it together with the data points.
• Calculate your mean square error.
• Now use the command polyfitcommand to find the best linear fit.
• Calculate with polyval to calculate the estimated data point and the error
• Plot the curve calculated by Matlab and its delta-values as error bars together with the data points in a figure.
• How good was your estimate?
• An alternative way to estimate the error is to use [p,S]=polyfit(x,y,n)returns the value S.normr. This value indicates the magnitude (Euclidean norm) of the residuals. (You can simply calculate it, for example, with residue_amount=S.normr for example).

*T8A5) The polyval function is not really needed for our application. Write a function with

• Input argument: Vector p of the coefficients and
• input argument: x-vector
• Return argument: the resulting vector of estimated values

• Using the last task as an example, test whether your estimates correspond to those of the polyval function.

T8A6) In order to better characterise the neuron, a series of measurements with current stimuli between 0 and 3.5nA was carried out and the response rate was measured 10 times in each case. The stimuli are stored under [stimT8.mat], the resulting rates under [rateT8.mat].

• Load the stimuli and rates and plot the rates (responses) as individual data points against the stimulus strengths.
• Find the best linear fit for the data and plot it.
• Extrapolation: What rate do you expect in response to a 5nA stimulus? Which one at -2nA?

T8A7) To check our extrapolation, the measurement series was extended The corresponding results are stored under [allereizeT8.mat] and [allatenT8.mat].

• Plot allrate against allereize as individual points together with the linear approximation already created for the part of the data processed in T8A3).
• Create a new linear approximation for all data points.
• Will the approximation be better if you use a second-order polynomial?
• What does it look like for polynomials of higher degree? To estimate the quality of the approximation, use
• • use both the mean square deviation of the measured points from the estimated curve
  • as well as S.normr
• Vary the degree of the polynomial systematically and record the two resulting errors in order to estimate how the fit improves as the number of coefficients increases.
• How many coefficients are needed to fit the data with polyfit satisfactorily?
• • When does the estimate stop improving significantly when additional coefficients are added?
  • Does the estimate get worse again at some point?
  • Do both error measures (mean square deviation and S.normr) lead to the same result?

T8A8) Calculate the mean and standard deviation for each of the current values with which measurements were made and plot this in a figure using errorbar. (Note: for this task it is more favourable to make a matrix from the long vector of measured values. This can be done with reshape).

*T8A9) Those interested in mathematics could take a look at the keyword "least squares fitting" to see how data fitting is carried out in Matlab and what other options the Curve Fitting Toolbox offers.

The aim of many experiments is to find a dependency between parameter values and a measured variable. An attempt is made to describe this dependency mathematically using a formula. This makes it possible to predict which measured values are expected for given parameter values. This procedure is also known as simulation. The experiment can then be used to check whether the predictions correspond to the actual measurements. If this is not the case, the model (i.e. the formula describing the dependency) must be adjusted accordingly.

Here we will use the sigmoid function as an example of a model. Sigmoid functions play an important role as saturation functions in biology, as there are hardly any examples of unbounded linear dependencies. Sigmoid functions can be represented using the formula
y=a./(1+(exp(b*(c-x))))
.

Tasks:

T8B1) In fact, the data just used are not real measurements, but points superimposed with noise on a sigmoid function:

• Create a vector xthat runs from -10 to 10 in small steps.
• Assign the parameters of the sigmoid formula with the values a=120; b=1; c=2. Calculate y according to the sigmoid formula.
• Plot the sigmoid function together with all rates and your best fit.
• How similar are the two lines and the measurement points?
• *) What standard deviation did I use to generate the data points?

T8B2) It is worth getting to know the sigmoid function in more detail. Write a programme:

• the three parameters of the sigmoid function should be systematically varied.
• Display the effect of each of the three parameters in a separate graph window as a set of curves.
• Consider: What effect do the three parameters have?
• (Tip: start from the values a=1, b=1, c=1 to see the effects clearly).

Biological processes are almost always sequences in time in which the current state depends on the past. In addition, chance normally plays a major role. Unfortunately, we do not have time to go into this topic in detail in this course. (The ICBM offers special courses on this topic.) Therefore, here is just an example of a simulation of such a time-dependent random process:

T8C1) Five rats live in a cage under permanent illumination. In our observation, we are currently only interested in the question of how many of the rats are asleep and how many are awake at a given time.

Write a simulation of the rats' sleep-wake behaviour in which you determine how many rats are awake every 10 minutes based on the following rules:

• Initial state: 3 rats are awake, 2 are asleep
• If only 0 or 1 rat is awake at a given time, the probability for each of the rats to be asleep at the next measurement is 80%.
• If 2 or more rats are awake, the probability for each rat is only 30% that it will be asleep at the next measurement (because the animals disturb each other).

• Run your simulation over a long period of time and record the number of rats awake.
• What is the percentage of time the animals spend asleep?

T8C2) Now consider rats in individual cages: 25 rats are housed in individual cages placed in five rows next to each other, resulting in a 5x5 matrix of cages.

• A random decision is made for each rat as to whether it is sleeping or not. The probability of sleeping is 40%.
• Simulate this random sleeping behaviour in time steps of 10 minutes for one week.
• Create a "film" in which the cages in which the rats are sleeping and those in which they are awake are shown graphically for each of your measurements (sleeping and awake are represented by two different colours).

T8C3) Now the rules become a little more interesting:

• Of the 25 rats, only those in neighbouring cages disturb each other while sleeping.
• If all neighbours are asleep, the probability for each rat is 80% that it will be asleep at the next measurement.
• If one neighbour is awake, the probability drops to 40%.
• If two neighbours are awake, the probability drops to 30%.
• If more than two neighbours are awake, there is only a 20% probability that the rat will sleep in the next time step.

• Run your simulation over a long period of time and watch the changing patterns as a "film".
• Record which of the rats is awake and when.
• What percentage of time do the individual animals spend asleep?
• Are there differences between the animals? Are these significant?

Command reference: Matlab syntax

T8H1) You inoculate 20 Petri dishes with 100 protozoa each and count the number of protozoa in the Petri dish 10 times per hour for 10 hours. This gives you the following matrix [populations.mat].

• Look at the measurement points graphically.
• Determine the polynomial that best describes the growth of the population.

*T8H2) In neuroscience, it is a common approach to consider the time of occurrence of individual action potentials as unimportant and only look at the spike rate (i.e. the number of action potentials occurring in a given time window). If the occurrence of individual spikes is then to be simulated, this is done with a random process, the so-called Poisson process. A Poisson process has the property that at each time step it is decided independently of all other time steps whether an event occurs or not. For example, if a neuron fires at 10 spikes/s and a sampling rate of 1000Hz is used, the probability of a spike occurring at each time step is 1%.

a) Write a function that receives a spike rate, the sampling rate and the length of the desired simulated response trace as input values and returns a vector with the simulated responses of the neuron.

b) Extend your function so that it does not receive a fixed value for the spike rate as input, but a vector that specifies the spike rate for each point in time of the response to be simulated (you then no longer need to pass the length, it results from the length of the vector of spike rates).

*T8H3) The Poisson process is the simplest conceivable neuron model and cannot explain many properties of biological nerve cells. A somewhat more realistic (but still very simple) model for action potential responses of a neuron is the integrate-and-fire neuron model. It is a differential equation model that describes the temporal change of the membrane voltage in combination with a threshold that determines the timing of the triggering of action potentials.

The demo programme [intfire.m] shows how an integrate-and-fire neuron model reacts to an external input current.

a) How does the model neuron react to other injected currents?

• Try making the pulse stronger or weaker.
• Make the pulse longer or shorter.
• What is the shortest pulse at which a spike is triggered?
• How does the model react to a sinusoidal stimulus?

b) What effect do the parameters have?

• What is the effect of the time constant?
• What effects does the firing threshold have?
• What happens if you start from a different membrane potential value?

c) How does noise affect the responses?

• Expand the model by adding the following to the integration at each time step tau_bin*randn during the integration. This step is intended to simulate the fact that nerve cells do not always respond in exactly the same way to the same stimulus.
• Run the simulation several times. How do the results change?
• Introduce a new parameter that scales the size of the noise and try out its effect.
• Modify the extended integrate-and-fire script to a function that returns the number of spikes generated but has no graphical output.
• Write a script that calls the new integrate-and-fire function 100 times and plots a histogram of the spike counts that occur.

**T8H4)Write a program that determines a sigmoid function for a given set of measurements that minimises the mean square distance to the measurements.

Test this function by creating artificial data in which a sigmoid function is superimposed with noise.

T8H5) Please review your notes from the course so that you can ask questions again tomorrow!

I would also like to ask you to read the additional texts"Concepts","Programming" and"Commands" and to think of questions for them as well. (Even if you may have done this before, they are probably easier to understand now).

If you have any comments on the script, the texts and tasks (what was incomprehensible, misleading, wrong or incomplete) I would be very grateful for a corresponding list!!!

And please remember to fill in the feedback form!