/*
Copyright (C) 2022 Leiden University Medical Center
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <APqr7.h>
#include <math.h>
#include <vector>
/*
*********
* APqr7 *
*********
This software provides an Action Potential Cure (APqr) to correct divergent
membrane potentials in excitable biological systems. The first X action
potentials are logged when the software starts, after which AP correction
starts from the (X+1)-st AP onwards. This correction occurs with the use of
dynamic pacth-clamp currrent injection.
IN:
*) Cm Capacitance of the cell
*) V_cutoff Threshold potential for the detection of the beginning
of an AP
*) Slope_tresh Slope threshold that defines the beginning of the
AP (mV/ms)
*) BCL_cutoff Threshold value for the end of an AP, given as a
percentage of the total APD
*) lognum Number of APs that need to be logged as a reference
*) Rm Initial resistance
*) Rm_corr_up Factor to increase Rm with when necessary
*) Rm_corr_down Factor to decrease Rm with when necessary
*) noise_tresh The noise level that is allowed around the ideal value
before correcting
OUT:
*) Vout voltage that is used to inject the calculated amount
of current into the excitable system
*/
/*
createRTXIPlugin
----------------
Creation of a new RTXI Plugin
IN:
*) None
OUT:
*) RTXIPlugin
*/
extern "C" Plugin::Object *createRTXIPlugin(void)
{
return new gAPqr7();
}
/*
vars[]
----
This is not a function, but rather the construction of a list. This list contains
all the variables that are visible and/or modifiable in the GUI of the software
module. There is the choice between INPUT, OUTPUT, PARAMETER and STATE.
INPUT: connected to the input port
OUTPUT: connected to the output port
PARAMETER: modifiable variable in the code
STATE: non-modifiable variable in the code
*/
static DefaultGUIModel::variable_t vars[] = {
{ "Vm (mV)", "Membrane potential (mV)", DefaultGUIModel::INPUT, },
{ "Iout (pA)", "Output current (pA)", DefaultGUIModel::OUTPUT, },
{ "Cm (pF)", "pF", DefaultGUIModel::PARAMETER
| DefaultGUIModel::DOUBLE, },
{ "V_cutoff (mV)", "Threshold potential for the detection of the beginning of an AP, together with Slope_thresh",
DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
{ "Slope_thresh (mV/ms)", "Slope threshold that defines the beginning of the AP (mV/ms)",
DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
{ "BCL_cutoff (pct)", "Threshold value for the end of an AP, given as a percentage of the total APD",
DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
{ "noise_tresh (mV)", "The noise level that is allowed before correcting", DefaultGUIModel::PARAMETER
| DefaultGUIModel::DOUBLE, },
{ "Rm (MOhm)", "MOhm", DefaultGUIModel::PARAMETER
| DefaultGUIModel::DOUBLE, },
{ "lognum", "Number of APs that need to be logged as a reference", DefaultGUIModel::PARAMETER
| DefaultGUIModel::DOUBLE, },
{ "Rm_corr_up", "To increase Rm when necessary", DefaultGUIModel::PARAMETER
| DefaultGUIModel::DOUBLE, },
{ "Rm_corr_down", "To decrease Rm when necessary", DefaultGUIModel::PARAMETER
| DefaultGUIModel::DOUBLE, },
{ "Correction (0 or 1)", "Switch Rm correction off (0) or on (1)",
DefaultGUIModel::PARAMETER | DefaultGUIModel::DOUBLE, },
{ "Period (ms)", "Period (ms)", DefaultGUIModel::STATE, }, // To check that the period taken by the algorithm is the same as the one i nthe control panel module
{ "Time (ms)", "Time (ms)", DefaultGUIModel::STATE, }, // To check that the algorithm is running
{ "APs2", "APs", DefaultGUIModel::STATE, }, // To check whether APs are being logged and the counter increases
{ "BCL2", "BCL", DefaultGUIModel::STATE, }, // To check what the eventual BCL of the ideal AP has become. You can see then if the APs were logged correctly
{ "act2", "0 or 1", DefaultGUIModel::STATE, }, // Switches from 0 to 1 and back continuously as a check to see whether you are computing error values and corrected values
};
/*
num_vars
--------
variable denoting the amount of variables that is displayed in the GUI
*/
static size_t num_vars = sizeof(vars) / sizeof(DefaultGUIModel::variable_t);
/*
gAPqr7
------
This function constructs the actual GUI by basing itself on the Default GUI Model.
It creates a module with a name, initializes the GUI, initializes the parameters,
adds a refresh, and allows you to resize.
IN:
*) None
OUT:
*) None
*/
gAPqr7::gAPqr7(void) : DefaultGUIModel("APqr7", ::vars, ::num_vars)
{
setWhatsThis(
"<p><b>APqr:</b><br>APqr7 </p>");
DefaultGUIModel::createGUI(vars, num_vars);
initParameters();
update(INIT);
refresh();
resizeMe();
}
gAPqr7::~gAPqr7(void) {}
/*
cleanup
-------
The APqr software makes use of three list structures which need cleaning after
a reset of parameters. The cleanup function takes care of this.
IN:
*) None
OUT:
*) None
*/
void gAPqr7::cleanup()
{
for(i=0;i<10000;i++){
Vm_log[i]=0;
Vm_diff_log[i]=0;
ideal_AP[i]=0;
}
}
/*
execute
-------
This is the main funtcion of the code that is looped through real-time.
It contains the four main parts of the algorithm:
1) Recording the ideal AP
2) Detecting AP upstrokes
3) Computing AP correction and outputting this
4) Updating the necessary variables
IN:
*) None
OUT:
*) Vout voltage that is used to inject the calculated amount
of current into the excitable system
*/
void gAPqr7::execute(void)
{
systime = count * period; // time in milli-seconds
Vm = input(0) * 1e2; // convert 10V to mV. Divided by 10 because
// the amplifier produces 10-fold amplified
// voltages. Multiplied by 1000 to convert
// V to mV.
Vm_log[count % (int)modulo] = Vm; // Logging the measured Vm in a list
// where the modulo component makes
// sure you keep cycling when you have
// reached the maximum number in the list.
// ****************************
// ****************************
// ** Recording the ideal AP **
// ****************************
// ****************************
if(count>(int)(1/period)-1 && (Vm - Vm_log[(count-(int)(1/period)) % (int)modulo]) >= slope_thresh && APs<lognum && enter == 0 && Vm > V_cutoff)
{
// This statement is entered whenever an upstroke is detected and the amount of
// recorded APs is smaller than lognum.
// The if conditions measure the following:
// 1) Whether you are far enough in the recording such that you don't accidentaly
// start in an ongoing AP
// 2) Whether two consecutive measuring points show a large enough slope that can
// be identified with an upstroke
// 3) Whether less than lognum APs were recorded
// 4) Whether you are currently not in an action potential
// 5) Whether the mesured voltage is above a voltage treshold
BCL = (APs==-1? 0: (BCL*APs + count2)/(APs+1)); // Rolling average of the basic cycle length
log_ideal_on = 1; // Switches on logging the AP
count2 = 0; // Resets the logging counter
enter = 1; // Switches on the indicator that an AP has started
APs++; // Counts the AP upstrokes that have passed
}
if((Vm - Vm_log[(count-(int)(1/period)) % (int)modulo]) < 0 && enter == 1)
{
// This statement is entered whenever the upstroke phase of an AP is over.
// The if conditions measure the following:
// 1) Whether two consecutive measuring points show a negative slope
// 2) Whether you currently are in an ongoing AP
enter = 0; // Switches off the indicator that an AP has started
}
if(APs<lognum && log_ideal_on == 1)
{
// This statement is entered whenever logging of the AP is on
// The if conditions measure the following:
// 1) Whether less than lognum APs were recorded
// 2) Whether the AP should be logged
ideal_AP[count2] = (ideal_AP[count2]*APs + Vm)/(APs+1); // Rolling average of the AP values
count2++; // Increasing the logging counter
}
// ****************************
// ****************************
// ** Detecting AP upstrokes **
// ****************************
// ****************************
if (act == 0 && (Vm - Vm_log[(count-(int)(1/period)) % (int)modulo]) >= slope_thresh && APs >= lognum && Vm > V_cutoff)
{
// This statement is entered whenever an upstroke is detected after the
// ideal APs have been recorded.
// The if conditions measure the following:
// 1) Whether currently nothing is being done or corrected
// 2) Whether two consecutive measuring points show a large enough slope that can
// be identified with an upstroke
// 3) Whether lognum APs were already recorded before
// 4) Whether the mesured voltage is above a voltage treshold
count = 0; // Reset the correction counter
act = 1; // Switch the correction on
}
// *************************************************
// *************************************************
// ** Computing AP correction and outputting this **
// *************************************************
// *************************************************
if (act == 1)
{
// This statement is entered whenever the instruction to correct the AP has
// been given.
Iout = Cm * (1/Rm) * (Vm - ideal_AP[count]); // Calculate the outward going current as
// a value proportional to capacitance,
// conductivity (1/resistance), and the error
output(0) = -Iout * 2.5e-3; // This is equal to Vout
// The factor 2.5e-3 comes from the conversion between current
// and voltage that is associated to the external command
// sensitivity of the Multiclamp 700B patch-clamp amplifier,
// which is 400 pA/V to be precise
Vm_diff_log[count] = Vm - ideal_AP[count]; // Log the errors
}
// **************************************
// **************************************
// ** Updating the necessary variables **
// **************************************
// **************************************
if(corr == 1 && act == 1 && count > 1 && abs(Vm_diff_log[count])>noise_tresh)
{
// This statement is entered whenever an update is needed in the resistance.
// The if conditions measure the following:
// 1) Whether correction adaptation is on
// 2) Whether currently there is correction going on
// 3) whether we are not in the very first step (gives errors)
// 4) whether the current error is larger than the noise threshold
if((Vm_diff_log[count-1] / Vm_diff_log[count]) < 0)
{
// This statement is entered whenever two consecutive error values have
// an opposite sign. This means that an overshoot in correction occurred.
// Therefore Iout should become less, and hence Rm should be increased.
Rm = Rm * Rm_corr_up; // Increase the resistance
}
if(abs(Vm_diff_log[count-1]) < abs(Vm_diff_log[count]) && (Vm_diff_log[count-1] / Vm_diff_log[count]) > 0)
{
// This statement is entered whenever two consecutive error values have
// the same sign, and when the error values increase in value. This means
// that the error is increasing. Hence we need to correct stronger and Iout
// should increase. As a consequence the resistance Rm should be decreased.
Rm = Rm / Rm_corr_down; // Decrease the resistance
}
}
if (count > BCL_cutoff*BCL)
{
// This statement is entered whenever the end of an AP is reached.
// The if condition measures the following:
// 1) Whether the current AP is further than a chosen cutoff of the pre-determined basic cycle length
act = 0; // Stop correcting during the last phase of the AP (is RMP)
output(0) = 0; // Send a 0 output since the last output is otherwise kept
}
count++; // End of the real-time loop, adjust the counter
}
/*
Update
------
This function updates the parameters of the code depending on the flag that is
given to it, where each flag is associated to a button.
INIT: associated to the loading of the module
MODIFY: associated to the Modify button
PERIOD: associated to the period linker with the "system control panel" module
PAUSE: associate to the pause button when pressing on it
UNPAUSE: associated to the pause button when unpressing it
IN:
*) flag Indicating the state of the update:
INIT, MODIFY, PERIOD, PAUSE, UNPAUSE
OUT:
*) None
*/
void gAPqr7::update(DefaultGUIModel::update_flags_t flag)
{
switch (flag)
{
case INIT:
setParameter("Cm (pF)", Cm);
setParameter("V_cutoff (mV)", V_cutoff);
setParameter("Rm (MOhm)", Rm);
setParameter("Rm_corr_up", Rm_corr_up);
setParameter("Rm_corr_down", Rm_corr_down);
setParameter("noise_tresh (mV)", noise_tresh);
setParameter("lognum", lognum);
setParameter("BCL_cutoff (pct)", BCL_cutoff);
setParameter("Slope_thresh (mV/ms)", slope_thresh);
setParameter("Correction (0 or 1)", corr);
setState("Time (ms)", systime);
setState("Period (ms)", period);
setState("APs2", APs);
setState("BCL2", BCL);
setState("act2", act);
break;
case MODIFY:
Cm = getParameter("Cm (pF)").toDouble();
Rm = getParameter("Rm (MOhm)").toDouble();
lognum = getParameter("lognum").toDouble();
V_cutoff = getParameter("V_cutoff (mV)").toDouble();
BCL_cutoff = getParameter("BCL_cutoff (pct)").toDouble();
noise_tresh = getParameter("noise_tresh (mV)").toDouble();
Rm_corr_up = getParameter("Rm_corr_up").toDouble();
Rm_corr_down = getParameter("Rm_corr_down").toDouble();
slope_thresh = getParameter("Slope_thresh (mV/ms)").toDouble();
corr = getParameter("Correction (0 or 1)").toDouble();
systime = 0;
count = 0;
APs = -1;
BCL = 0;
log_ideal_on = 0;
enter = 0;
count2 = 0;
cleanup();
break;
case PERIOD:
period = RT::System::getInstance()->getPeriod() * 1e-6; // time in milli-seconds
modulo = (1.0/(RT::System::getInstance()->getPeriod() * 1e-6)) * 1000.0;
break;
case PAUSE:
output(0) = 0.0;
Iout = 0;
act = 0;
systime = 0;
break;
case UNPAUSE:
break;
default:
break;
}
}
/*
initParameters
--------------
This function sets all values to their defaults when no external parameters are provided
through the GUI interface.
IN:
*) None
OUT:
*) None
*/
void gAPqr7::initParameters()
{
Vm = -80; // mV
Cm = 150; // pF
Rm = 150; // MOhm
slope_thresh = 5.0; // mV
corr = 1;
Iout = 0; // pA
output(0) = -Iout * 0.5e-3;
period = RT::System::getInstance()->getPeriod() * 1e-6; // ms
systime = 0;
count = 0;
act = 0;
Rm_corr_up=8;
Rm_corr_down=2;
noise_tresh = 0.5; // mV
BCL = 0; // ms
count2 = 0;
APs = -1;
V_cutoff = -40; // mV
BCL_cutoff = 0.98;
enter = 0;
log_ideal_on = 0;
lognum = 3;
modulo = (1.0/(RT::System::getInstance()->getPeriod() * 1e-6)) * 1000.0;
}