QUALITY CONTROL OF MEDICINES

Development and validation of a method for authentication of a encapsulated medicinal product by the non-destructive express method of NIR spectrometry and chemometrics analysis

Author information

1 — The Federal State Budget Institution “Information center for expertise, accounting and analysis of circulation of medical products” of Federal Service for Surveillance in Healthcare, 4, bld.1, Slavyanskaya Square, 109012, Moscow, Russian Federation.

ORCID: https://orcid.org/0000-0002-2829-8801

2 — The Federal State Budget Institution “Information center for expertise, accounting and analysis of circulation of medical products” of Federal Service for Surveillance in Healthcare, 4, bld.1, Slavyanskaya Square, 109012, Moscow, Russian Federation.

3 — The Federal State Budget Institution “Information center for expertise, accounting and analysis of circulation of medical products” of Federal Service for Surveillance in Healthcare, 4, bld.1, Slavyanskaya Square, 109012, Moscow, Russian Federation.; Semenov Federal Research Center for Chemical Physics RAS, 4, Kosygin str. 119334 Moscow, Russian Federation.

ORCID: https://orcid.org/0000-0002-0146-8284

4 — The Federal State Budget Institution “Information center for expertise, accounting and analysis of circulation of medical products” of Federal Service for Surveillance in Healthcare, 4, bld.1, Slavyanskaya Square, 109012, Moscow, Russian Federation.; Federal State Autonomous Educational Institution of Higher Education I.M. Sechenov First Moscow State Medical University of the Ministry of Health of the Russian Federation, 8 bld. 2, Trubetskaya St., Moscow 119048, Russian Federation.

ORCID: https://orcid.org/0000-0002-0636-3254

Published: 26.06.2024

One of the ways to protect the market from substandard and falsified medicines is the authentication of medicines. The combination of near-infrared spectrometry (NIR spectrometry) with chemometric analysis (NIR analysis) is one of the priority methods for this purpose, since it allows for rapid measurements without violating the integrity of the primary packaging and does not require sample preparation. The medicine “Ribavirin, capsules 200 mg” has been choosing as an object of the study. The goal of the work is development of a NIR spectral model for the authentication of the above medicine. The method of one-class classification has been using for medicines modeling. Various ways of preprocessing spectra and methods of calculating the model was studied. Validation of the models showed that the optimal condition for developing the model is the method of calculating the Pearson correlation coefficient without preprocessing the spectra in the operating range of 12 000–6000 cm-1 and 5500–4400 cm-1. The sensitivity of the models is 100.0 %, and the specificity is 99.9 %. The significance of the choice of the method of constructing the model, the working range of measurement and the composition of the medicine on the model specificity is shown. The analysis of the results showed that all factors (the operating range, the method of preprocessing the spectra, the method of developing the model and the spectral characteristics of the sample) are important significant in the process of developing the model and in different combinations can manifest themselves with different strengths. Therefore, the authentic model can be developed only by studying each options.

Keywords: medicines authentication, NIR spectroscopy, chemometric analysis, method of one-class classification, spectra preprocessing, validation of the medicines mathematic model

Background.

Одним из способов защиты рынка от некачественных и фальсифицированных лекарственных препаратов является их аутентификация. В литературе описаны разные способы ее реализации [1–7], среди них наиболее часто используется спектрометрия в ближней инфракрасной области (БИК-спектрометрия) в сочетании с хемометрическими методами анализа (БИК-анализ).

Очевидными достоинствами БИК-спектрометрии являются возможность выполнения испытания препарата без разрушения первичной упаковки и самого препарата, быстрота измерения, отсутствие деструктивного воздействия на анализируемый образец, что позволяет после завершения испытания вернуть его в гражданский оборот для дальнейшей реализации.

В отличие от идентификации, целью которой является подтверждение подлинности действующего вещества в препарате, аутентификация методом БИК-спектрометрии нацелена на создание спектрального образа препарата (БИК-спектральная модель), в котором участвуют как действующие, так и вспомогательные вещества. Поэтому важным этапом аутентификации лекарственных средств БИК-анализом является построение БИК-спектральной модели препарата, которая бы позволила отличить его от «чужих образцов», а также выявить изменения, произошедшие с ними либо в процессе производства, либо хранения.

Цель работы

Целью настоящей работы является разработка математической модели для аутентификации капсулированного лекарственного препарата неразрушающим экспресс-методом БИК-спектрометрии и хемометрическими методами анализа для своевременного выявления на рынке некачественных и фальсифицированных образцов препарата.

Объекты и методы исследования

В качестве объекта исследования были выбраны капсулы препарата «Рибавирин» дозировкой 200 мг производителя Пс. Предоставленные производителем образцы были упакованы в контурную ячейковую упаковку из поливинилхлорида (ПВХ) и алюминиевой фольги (первичная упаковка) и картонную пачку (вторичная упаковка).

Препарат «Рибавирин» широко используется в качестве противовирусного средства в основном для лечения гепатита С и вирусных геморрагических лихорадок, часто назначается для лечения орфанных заболеваний.

Спектры в ближней ИК области (БИК-спектры) были получены на ИК Фурье спектрометре MPA фирмы «Брукер» (Германия) с помощью оптоволоконного датчика методом диффузного отражения. Разрешение – 8 см-1, количество сканов – 16, область измерения – от 12 000 см-1 до 4000 см-1, фазовое разрешение – 32, интерполяция – 2, базовая линия – по эталону из тефлона.

Построение и валидацию модели осуществляли с помощью подпрограммы «Ident» из программы «Opus 8.5» фирмы «Брукер» (Германия).

Учитывая, что работа была выполнена на разных единицах указанного прибора, для исключения различия в шкале волновых чисел используемых приборов все БИК-спектры были совмещены с условно стандартным спектром.

В качестве математического метода построения модели был выбран метод одноклассовой классификации, основной задачей которого является формирование целевого класса в соответствии со свойствами исследуемой группы объектов [1–2].

В соответствии с требованиями ОФС.1.2.1.1.0001.15 «Спектрометрия в ближней инфракрасной области»1 модель, построенную методом одноклассовой классификации, валидировали на чувствительность (степень принятия образцов целевого класса) и специфичность (степень отклонения «чужих образцов», т.е. не принадлежащих к целевому классу).

Чувствительность модели оценивали в два этапа: с помощью спектров, используемых для построения модели (внутренняя чувствительность), и спектров, полученных с дополнительных образцов того же препарата, но неиспользованных в построении модели (внешняя чувствительность).

Оценку модели на специфичность провели на спектрах поверочных образцов, в качестве которых были использованы выпущенные в гражданский оборот образцы капсул «Рибавирин» с той же дозировкой и в такой же упаковке, как исследуемый образец, но других производителей (табл. 1, рис. 1). В работе также были использованы спектры образцов препарата производителя Пк нового состава (Пк-н) и предыдущего (Пк-с).

Из каждой серии предоставленных производителями образцов исследуемого препарата методом случайной выборки были отобраны по 12 капсул: 10 – для построения модели и валидации на внутреннюю чувствительность, 2 – для ее валидации на внешнюю чувствительность. Для каждой капсулы были получены по 3 БИК-спектра с плотным прижатием датчика к капсуле на уровне ее корпуса через прозрачную часть первичной упаковки.

Дополнительно были получены спектры капсул с содержимым без упаковки, самой капсулы и упаковки (рис. 3). Для этого из набора для построения модели из каждой серии были выбраны по 2 капсулы, удалены из упаковки и получены по 3 спектра на уровне корпуса капсулы.

БИК-спектры пустых капсул (заполненных воздухом) получались нестабильными. Для устранения этого эффекта был подобран спектрально инертный наполнитель – субстанция натрия хлорид с массовой долей действующего вещества не менее 99,9 %. Однако в БИК-спектре указанной субстанции были обнаружены две широкие полосы, совпадающие с полосами воды (рис. 1). Для удаления абсорбционной воды образец натрия хлорида высушили при 105 °С в течение 1 часа: полосы в спектре не исчезли (рис. 1). Только прокаливанием субстанции при 650 °С в течение 1 часа удалось избавиться от указанных полос в спектре (рис. 1). Это свидетельствует о том, что в субстанции натрия хлорида присутствовала кристаллизационная вода.

Спектры пленки из поливинилхлорида (ПВХ), входящей в состав блистера, были получены в трех повторах для каждой упаковки из каждой серии препарата путем плотного прижатия датчика к поверхности полимерной части упаковки между ячейками.

Необходимые для анализа спектра испытуемого препарата БИК-спектры входящих в него действующего и вспомогательных веществ (рибавирин производства ЗАО «Активный компонент»; лактозы моногидрат производства компании «DFE Pharma», Нидерланды; магния стеарата фирмы «Tablube Nitica», Индия; кремния диоксида коллоидальный компании «Evonik Industries AG», Германия) были также получены в трех повторах путем прямого измерения субстанций во флаконах.

Ключевой позицией в моделировании является построение и валидация математической модели, максимально адекватной составу образцов лекарственного препарата, учитывающей особенности технологии производства и отличающей «чужие образцы». Для решения поставленной цели, с учетом имеющегося оборудования и программного обеспечения, оптимальными являются классификационные методы математического моделирования и одноклассовые по способу их построения [1–2].

Для построения одноклассовой модели в программе OPUS предложены два метода вычисления: расстояние Евклида (стандартный метод) и коэффициент корреляции Пирсона. В качестве предобработки спектров предусмотрены векторная нормализация, первая или вторая производные с шириной окна от 5 до 25 точек или сочетание производных с векторной нормализацией.

В том и другом случае пороговое значение модели (DT) представляет собой сумму максимального расстояния (Dmax) между средним спектром и спектрами образцов, используемых для построения модели (ядро модели), и произведение стандартного отклонения (S0) и коэффициента (Кх):

DT = Dmax + КхS0.

По умолчанию коэффициент Кх имеет значение 0,25, но его можно изменить от 0 до 1, передвигая таким образом границу модели с целью ее оптимизации.

Специфика коэффициента корреляции Пирсона, в отличие от расстояния Евклида, состоит в том, что в процессе вычисления параметров модели одновременно выполняется нормализация спектров, поэтому при его использовании в качестве метода вычисления модели дополнительная предобработка спектров в виде векторной нормализации не проводится.

Существенным фактором для построения модели является выбор рабочего диапазона, который на первом этапе проводится визуально путем анализа и выявления тех областей спектра, в которых проявляются полосы действующего и вспомогательных веществ.

Визуальный анализ среднего БИК-спектра образца препарата показал, что в нем превалируют полосы желатина (рис. 1), поскольку его содержание в препарате составляет около 30 % и он локализован в верхнем слое препарата.

Действующее вещество (рибавирин) проявляется малоинтенсивными полосами при волновых числах 6766, 6369, 5052, 4731, 4612, 4528 и 4482 см-1. Лактоза, которая является основным наполнителем в содержимом капсулы (32 %), проявляется в виде широкой полосы в области волновых чисел 6600–6488 см-1, плеча в области 5183–5142 см-1 и в виде полосы при волновом числе 4015 см-1 (рис. 1). Кремния диоксид не имеет характерных полос в БИК-спектре, а полос магния стеарата при содержании его в препарате в количестве менее 1% практически не видно. ПВХ блистера проявляется в виде трех интенсивных полос при волновых числах 5828, 4330, 4215 см-1, малоинтенсивной полосы 5732 см-1 и плеча 5702 см-1.

С точки зрения расположения полос компонентов препарата область от 12 000 см-1 до 9000 см-1 является малоинформативной, в ней часто проявляются спектральные шумы, поэтому она могла бы быть удалена из исследования. Но известно, что математические методы могут извлекать из спектра невидимую глазом информацию.

В области 4150–4000 см-1 проявляются искажения спектров за счет нестабильной работы оптоволоконного датчика [1], а в областях от 5500–6000 см-1 и от 4400– 4050 см-1 оказывает влияние материал упаковки.

Присутствующая в образцах вода часто приводит к ухудшению разделяющих способностей моделей за счет увеличения максимального расстояния между спектрами и стандартного отклонения, поэтому область основного пика воды (5300– 4950 см-1) желательно было бы удалить. Однако в данном случае это является нецелесообразным, так как в этой области присутствуют полосы действующего и вспомогательных веществ препарата.

Таким образом, для построения моделей были выбраны четыре рабочих спектральных диапазона:

1) 12 000–4400 см-1;

2) 9000–4400 см-1;

3) 12 000–6000+5500–4400 см-1;

4) 9000–6000+5500–4400 см-1.

Перед построением моделей необходимо было убедиться в том, что испытуемые образцы принадлежат к одной выборке и в них отсутствуют грубые выбросы. Для этой цели была использована перекрестная валидация, суть которой состоит в исключении спектров образцов одной серии из модели, построении новой модели на спектрах образцов оставшихся серий и в оценке степени приятия этой моделью исключенных образцов. Исследование однородности выборки проводили в наиболее широком рабочем диапазоне (12 000–4400 см-1), в качестве предобработки спектров использовали векторную нормализацию, для построения модели – стандартный метод вычисления расстояния Евклида: все спектры были приняты моделями, что свидетельствует об однородности выборки.

Для выбора оптимальных условий построения БИК-спектральной модели были проанализированы следующие сочетания способов вычисления модели и предобработки спектров:

1) расстояние Евклида, без предобработки (РЕ бп);

2) расстояние Евклида, векторная нормализация (РЕ вн);

3) расстояние Евклида, первая производная со скользящим окном длиной 5 точек сглаживания (РЕ 1пр5);

4) расстояние Евклида, первая производная со скользящим окном длиной 25 точек сглаживания (РЕ 1пр25);

5) расстояние Евклида, вторая производная со скользящим окном длиной 5 точек сглаживания (РЕ 1пр5);

6) расстояние Евклида, вторая производная со скользящим окном длиной 25 точек сглаживания (РЕ 1пр25);

7) расстояние Евклида, первая производная со скользящим окном длиной 5 точек сглаживания плюс векторная нормализация (РЕ 1пр5+вн);

8) расстояние Евклида, первая производная со скользящим окном длиной 25 точек сглаживания плюс векторная нормализация (РЕ 1пр25+вн);

9) расстояние Евклида, вторая производная со скользящим окном длиной 5 точек сглаживания плюс векторная нормализация (РЕ 1пр5+вн);

10) расстояние Евклида, вторая производная со скользящим окном длиной 25 точек сглаживания плюс векторная нормализация (РЕ 1пр25+вн);

11) коэффициент корреляции Пирсона, без предобработки (КК бп);

12) коэффициент корреляции Пирсона, первая производная со скользящим окном длиной 5 точек сглаживания (КК 1пр5);

13) коэффициент корреляции Пирсона, первая производная со скользящим окном длиной 25 точек сглаживания (КК 1пр25);

14) коэффициент корреляции Пирсона, вторая производная со скользящим окном длиной 5 точек сглаживания (КК 2пр5);

15) коэффициент корреляции Пирсона, вторая производная со скользящим окном длиной 25 точек сглаживания (КК 2пр25).

В основе сравнения испытуемого образца и модели лежит стандартный метод вычисления спектрального расстояния (D) между полученным спектром со средним спектром модели путем вычисления евклидового расстояния по всем точкам спектрального диапазона модели [1].

В случае аутентичности исследуемого образца его спектральное расстояние (D) не должно превышать порогового значения модели (DT).

Валидация построенных моделей на внутреннюю и внешнюю чувствительность показала, что модели полностью принимают спектры своих образцов (табл. 2).

Результаты валидации модели на специфичность представлены в таблице 2, где для каждого производителя указано количество спектров, отвергнутых моделью и принятых ею.

Вычисление общей специфичности модели проводили путем суммирования отвергнутых моделью спектров по всем проверочным образцам, делением полученной суммы на общее количество спектров и умножением на 100 для перевода полученного результата в процентное соотношение (табл. 2, рис. 4). Аналогичным образом была вычислена общая чувствительность модели по принятию своих образцов, которая во всех случаях была равна 100,0 % (табл. 2).

Для определения степени влияния образцов препарата каждого производителя на общую специфичность были построены диаграммы специфичности модели по производителям (рис. 5–10). С целью исключения влияния количества используемых в анализе спектров на специфичность (табл. 2) было вычислено процентное отношение выбросов к количеству используемых по каждому производителю спектров.

Из гистограммы общей специфичности (рис. 4) видно, что наилучшие результаты достигнуты при использовании рабочего диапазона измерения 12 000–6000+5500–4400 см-1, вычисление модели с помощью коэффициента корреляции Пирсона без предобработки: общая специфичность модели составила 99,9%.

Немного хуже был получен результат при вычислении модели с помощью расстояния Евклида с предварительной предобработкой спектров методом векторной нормализации: 99,5%.

Аналогичные результаты по оптимальным условиям видны на диаграммах специфичности моделей по производителям (рис. 5–10).

Увеличение специфичности модели за счет уменьшения ее порогового значения DT с помощью коэффициента Кх оказалось нецелесообразным, так как величина их качества совпадений была меньше максимального расстояния между спектрами Dmax, т.е., по сути, выбросы расположились в ядре модели.

Анализ гистрограммы общей специфичности моделей выявил следующие закономерности:

  • модели, построенные методом расстояния Евклида и с использованием векторной нормализации в качестве предобработки спектров, имеют более высокую специфичность, чем модели без предобработки;
  • используемые в качестве предобработки производные нивелируют различия между спектрами и, как следствие, уменьшают специфичность модели при обоих способах ее построения;
  • увеличение в производных длины скользящего окна с 5 точек сглаживания до 25 увеличивает специфичность;
  • сочетание производных и векторной нормализации уменьшает специфичность модели при обоих способах ее построения;
  • широкий рабочий диапазон показывает лучшие результаты по специфичности модели, чем тот диапазон, в котором сосредоточены все полосы компонентов препарата, за исключением использования при предобработке спектров второй производной с длиной скользящего окна 5 точек;
  • удаление области расположения полосы ПВХ блистера из рабочего диапазона измерения приводит к увеличению специфичности модели только в случае ее построения без предобработки спектров или с использованием векторной нормализации; при использовании предобработки спектров с помощью производных специфичность модели снижается.

Диаграммы специфичности моделей по производителям (рис. 5–10) показали, что наиболее близким по спектру к анализируемому образцу препарата является образец препарата Пк-н, спектры которого часто полностью принимаются моделью (рис. 7). Образцы препарата По и Пф в большинстве случаев не принимаются моделями, т.е. имеют значительные спектральные отличия от образцов препарата Пс, на которых были построены модели (рис. 9 и 10).

Анализ диаграмм специфичности моделей по производителям показал следующее:

  • применение векторной нормализации в качестве предобработки спектров чаще приводит к увеличению специфичности модели, но может не повлиять или даже уменьшить ее, как это видно на рисунках 6, 8 и 9;
  • сочетание производных и векторной нормализации в большинстве случаев уменьшает специфичность модели или оставляет ее без изменения; в некоторых случаях видна зависимость этого сочетания от рабочего диапазона измерения (рис. 8–10);
  • использование производных для предобработки спектров приводит как к увеличению специфичности модели (рис. 5), так и к ее уменьшению (рис. 7 и 9);
  • для образца Пк-с изменение специфичности моделей при использовании производных зависит еще от диапазона измерения (рис. 8);
  • увеличение длины скользящего окна с 5 до 25 точек сглаживания при использовании производных в качестве предобработки спектров в большинстве случаев увеличивает специфичность, особенно эта закономерность четко прослеживается при вычислении модели коэффициентом корреляции Пирсона; изменение специфичности в ряде случаев зависит от рабочего диапазона измерения (рис.8);
  • замена метода вычисления модели с расстояния Евклида на коэффициент корреляции Пирсона приводит как к увеличению специфичности модели, так и к ее снижению, в большинстве случаев это зависит от диапазона измерения и состава поверочного образца (рис. 5–10);
  • наиболее индифферентным к диапазону измерения является образец препарата Пвф, на его диаграмме специфичности расхождения между диапазонами измерения наблюдаются только при применении второй производной в качестве предобработки спектров (рис. 6). В случае использования других образцов выбор оптимального диапазона измерения неоднозначен. Например, модель, построенная вычислением коэффициента корреляции Пирсона без предобработки в широких диапазонах измерения, отвергает спектры образцов препарата Пк-с, а в узких – принимает их. И наоборот, модель, построенная вычислением расстояния Евклида с использованием второй производной с длиной скользящего окна с пятью точками сглаживания в качестве предобработки, отвергает те же образцы в узком диапазоне измерения и принимает в широком (рис. 8). Такая же противоречивая картина наблюдается с использованием диапазонов измерения с удалением области полосы ПВХ упаковки и без нее (рис. 8–10).

В результате проведенного исследования общих для всех случаев закономерностей, однозначно увеличивающих специфичность модели, не выявлено. Это свидетельствует о том, что все факторы, а именно: рабочий диапазон, способ предобработки спектров, метод построения модели и спектральная характеристика образца, являются значимыми в этом процессе, и в разных сочетаниях могут проявляться с разной силой.

Conclusion.

На примере капсул «Рибавирин» показана процедура построения и валидации модели для аутентификации лекарственных препаратов в капсулированной форме методом БИК-спектрометрии в сочетании с хемометрическими методами анализа. Валидация моделей позволила определить оптимальные условия построения для образцов препарата «Рибавирин», капсулы, 200 мг производства Пс: рабочий диапазон измерения 12 000–6000+5500–4400 см-1, вычисление модели с помощью коэффициента корреляции Пирсона без предобработки. Чувствительность модели – 100,0%, общая специфичность –99,9%. Анализ гистрограммы общей специфичности моделей и диаграмм специфичностей моделей по производителям показал, что все факторы (рабочий диапазон, способ предобработки спектров, метод построения модели и спектральная характеристика образца) являются значимыми в этом процессе, и в разных сочетаниях могут проявляться с разной силой. Это является причиной возникновения разнообразных, часто противоположных результатов. В этих условиях поиск оптимальных условий построения аутентичной модели может быть осуществлен путем исследования всех возможных вариантов, как это описано в статье.

___________________________________________________________________

1 Статья Государственной фармакопеи Российской Федерации XIV издания.

  1. Titova A.V. , Rodionova O.Ye., Godin F.Yu., Balyklova K.S. Development and validation of a method for authenticating desloratadine tablets by NIR spectrometry // Vestnik Roszdravnadzora. – 2022. – Vol. 4. – P. 74–80. (in Russian).
  2. Rodionova, O.Ye. Qualitative and quantitative analysis of counterfeit fluconazole capsules: A non-invasive approach using NIR spectroscopy and chemometrics [Electronic Resource] / O.Y. Rodionova, A.V. Titova, N.A. Demkin, K.S. Balyklova, A.L. Pomerantsev // Talanta. – 2019; 195: 662–667. DOI: 10.1016/j.talanta.2018.11.088. – Access mode: https://pubmed.ncbi.nlm.nih.gov/30625598/ (date of request: 10.08.2023).
  3. Naughton B. Effectiveness of medicines authentication technology to detect counterfeit, recalled and expired medicines: a two-stage quantitative secondary care study [Electronic Resource] / B. Naughton, L. Roberts, S. Dopson, S. Chapman, D. Brindley // BMJ Open. – 2016; 6(12): e013837. DOI: 10.1136/bmjopen-2016-013837. – Access mode: https://pubmed.ncbi.nlm.nih.gov/27940634/ (date of request: 10.08.2023).
  4. Assi S. Authentication of Antibiotics Using Portable Near- Infrared Spectroscopy and Multivariate Data Analysis [Electronic Resource] / S. Assi, B. Arafat, Kathryn Lawson-Wood, I. Robertson // Applied Spectroscopy. – 2021; 75(4): 434– 444. DOI: 10.1177/0003702820958081. – Access mode: https://www.researchgate.net/publication/343845424_EXPRESS_Authentication_of_Antibiotics_Using_ Portable_Near-Infrared_Spectroscopy_and_Multivariate_Data_Analysis (date of request: 10.08.2023).
  5. Chen H., Lin Z., Tan C. Application of near-infrared spectroscopy and class-modeling to antibiotic authentication / H. Chen, Z. Lin, C. Tan // Analytical Biochemistry. – 2020; 590:113514. DOI: 10.1016/j.ab.2019.113514. (date of request: 10.08.2023).
  6. Chen H. Express detection of expired drugs based on nearinfrared spectroscopy and chemometrics: A feasibility study [Electronic Resource] / H. Chen, C. Tan, Z. Lin // Spectrochimica acta. Part A: Molecular and biomolecular spectroscopy. – 2019; 220: 117153. DOI: 10.1016/j.saa.2019.117153. – Access mode: https://pubmed.ncbi.nlm.nih.gov/31141774/ (date of request: 10.08.2023).
  7. Hattori Y. Rapid identification of oral solid dosage forms of counterfeit pharmaceuticals by discrimination using nearinfrared spectroscopy [Electronic Resource] / Y. Hattori, Y. Seko, J. Peerapattana, K. Otsuka, T. Sakamoto, M. Otsuka // Bio-Medical Materials and Engineering. – 2018; 29(1): 1–14. DOI: 10.3233/ BME-171708. – Access mode: https://pubmed.ncbi.nlm.nih. gov/29254069/ (date of request: 10.08.2023).