APqr7.cpp

/*
 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;
}