Literature DB >> 27224203

Femtosecond-Laser-Pulse Characterization and Optimization for CARS Microscopy.

Vincenzo Piazza1, Giuseppe de Vito1,2, Elmira Farrokhtakin3, Gianni Ciofani3,4, Virgilio Mattoli3.   

Abstract

We present a simple method and its experimental implementation to determine the pulse durations and linear chirps of the pump-and-probe pulse and the Stokes pulse in a coherent anti-Stokes Raman scattering microscope at sample level without additional autocorrelators. Our approach exploits the delay line, ubiquitous in such microscopes, to perform a convolution of the pump-and-probe and Stokes pulses as a function of their relative delay and it is based on the detection of the photons emitted from an appropriate non-linear sample. The analysis of the non-resonant four-wave-mixing and sum-frequency-generation signals allows for the direct retrieval of the pulse duration on the sample and the linear chirp of each pulse. This knowledge is crucial in maximizing the spectral-resolution and contrast in CARS imaging.

Entities:  

Mesh:

Year:  2016        PMID: 27224203      PMCID: PMC4880195          DOI: 10.1371/journal.pone.0156371

Source DB:  PubMed          Journal:  PLoS One        ISSN: 1932-6203            Impact factor:   3.240


Introduction

In the coherent-anti-Stokes-Raman-scattering (CARS) process a pair of photons (“pump” and “Stokes”) is exploited to excite selected vibrational levels of a sample by choosing them so that their frequency difference matches the vibrational frequencies of interest. A third photon, from the same source as the pump in the so-called two-color CARS implementations, coherently probes the phonon population of the vibrational mode generating strong anti-Stokes emission [1]. This approach is being successfully used to obtain video-rate label-free imaging of biological samples with chemical selectivity [2] and to probe the orientation of chemical bonds in biological samples [3-6]. The third-order nature of the process suggests using ultrashort laser pulses to maximize the generated signal intensity at constant average power delivered to the sample. Conversely, reducing the pulse durations below the lifetime of the Raman modes under investigation (~1 ps in typical biological samples) yields reduced contrast due to non-resonant four-wave-mixing generation as well as poor spectral selectivity. To overcome this issue, ps laser sources are routinely used. On the other hand, these sources are not particularly efficient in multimodal microscopes where–apart from CARS signals–one is interested in generating and detecting second or third harmonic or two-photon fluorescence [5-11], where pulse durations in the 100-fs range are desirable. A flexible approach consists in introducing a controlled linear chirp in both beams and increase the duration of the pulses. This can be achieved using pulse stretchers or dispersing elements (e. g. blocks of glass) [12, 13] placed in the paths of ~100-fs pump and Stokes beams. If the chirp is chosen appropriately, the instantaneous frequency difference between the pump and Stokes beam is constant as a function of time. This approach was shown to yield improved spectral selectivity, comparable with that obtainable with transform-limited pulses with the same duration of the chirped pulses, and increased ratio between the resonant and non-resonant signals. It was also exploited in stimulated-Raman scattering [14] and implemented using blocks of glass as stretchers [15] also in single-source CARS microscopy [16], and in broadband CARS [17]. Removing the chirping elements easily restores the sub-ps duration of the pulses if needed. A similarly flexible approach was also proposed in Ref. [18] where only the pump-and-probe beam was strongly chirped. A successful implementation of these techniques requires determining the characteristics of the pulses at the sample level in order to appropriately match their chirps and durations. Pulses are routinely characterized by means of autocorrelators, where a time-delayed fraction of a pulse interferes on an object having a non-linear optical response with the non-delayed remaining of the same pulse. Second-order processes are recorded as a function of the time delay yielding the beam autocorrelation function and–when the spectral bandwidth and pulse-envelope shape are known–the pulse duration and chirp. Anyway, this characterization cannot be easily performed at the focal plane of the microscope objective and the additional contributions to the total group delay dispersion (GDD) are typically only estimated. The approach proposed here leads to distinctive advantages even when the contribution of the objective to the total GDD is small with respect to the other optical components in the microscope (or if it can be estimated). In this case the pulse characteristics could be determined with an external autocorrelator by removing the objective and the condenser but this process is cumbersome and can lead to misalignments in the optical path. It should also be noted that optimization strategies common for example in two-photon-excited-fluorescence microscopy, based on maximizing the signal by tuning a pre-chirp device, cannot be applied here. Indeed, the setting that would lead to maximum signal would likely be different to that leading to optimal spectral resolution due to the presence of non-resonant background. More sophisticated and precise laser-pulse-characterization techniques exist [19], e. g. based on frequency-resolved optical gating (FROG) [20] and improvements thereof. FROG is able to fully determine the intensity profile and phase of an unknown arbitrary pulse by optical gating a portion of a pulse with a delayed portion of the same while measuring the emitted-signal spectrum and using an iterative algorithm to analyze the collected spectra [21]. The signal may come from different types of non-linear interactions, e.g. second-harmonic generation, third-harmonic generation, self-diffraction, etc. [22] It is worth mentioning that the issue here of interest, i.e. the characterization of two unknown laser pulses–rather than just one,–was also extensively addressed in the past years [23] with specific variations of the FROG scheme, in particular with the so-called Blind FROG [24] and the more robust Double Blind FROG (DBFROG) [25]. In DBFROG, each of the two pulses is used as an optical gate for the other and two different spectra of the chosen non-linear signal (e.g. polarization gate) are then recorded for each position of the delay line. An iterative algorithm yields the full characterization of the two pulses. Of note, the typical CARS microscope with its delay line might be easily modified to realize a DBFROG-like geometry with the addition of an external spectrometer. Anyway, while the DBFROG technique and similar approaches are extremely powerful in the determination of unknown arbitrary pulses, they require complex iterative pulse retrieval algorithms. In the case of CARS microscopes, where one starts with known pulses and is interested in determining the additional linear chirp introduced to achieve spectral focusing and the total pulse duration at the sample, a simpler method would be desirable. In fact, simplified implementations are often used: in Ref. [26] XFROG SFG traces are recorded using a Fourier-transform-limited pump-and-probe pulse to characterize the Stokes pulse and vice-versa. In the Supplementary Information (S3 Fig) we shall exploit this approach to test the assumption of Gaussian beam profiles. Here a straightforward method to determine the main features of both the pump-and-probe pulse and the Stokes pulse at the sample is proposed and experimentally verified. By measuring the intensities of the emitted non-resonant four-wave-mixing (FWM) and sum-frequency-generation (SFG) signals as a function of their relative delay δt, the pulse durations at the sample level can be straightforwardly measured. Differently from conventional FROG-based approaches, this step does not even require to frequency resolve the emitted signal signals. Measuring the signal intensity with a PMT (a pair of filters to isolate the FWM and SFG will be needed) as a function of δt is enough to determine both pulse durations. If the pulse characteristics at the laser output are known, then it is also easy to calculate the GDDs at sample level. If the emitted signals can be spectrally analyzed, the GDD of each pulse can also be determined experimentally. The only assumptions behind this approach are that: i) the non-linear responses of the chosen medium are only weakly dependent on the wavelength of the excitation photons, at least within the bandwidths of the pulses, which can be easily verified experimentally. ii) Pump and Stokes beam powers are chosen in order to avoid pump depletion. Also this can be checked for experimentally (see for example S2 Fig). iii) In the following we shall assume that the pulses have a Gaussian profile. Interestingly, the approach detailed here does not require that one of the pulses is much longer than the other as in Ref. [26]. In fact, by collecting both the FWM and SFG signals also pulses of similar durations and GDDs can be characterized. It should be noted here that for the SFG signal to be present, a non-centrosymmetric sample is needed.

Theory

The electric field of a Fourier-limited Gaussian pulse having duration τ, amplitude A1/2, and centered around an angular frequency ω can be described as: After travelling through a dispersive medium the pulse acquires a GDD, quantified here by a parameter a. Neglecting third-order dispersion (TOD) terms and higher (an estimation of TOD terms is presented in S1 Text), the new pulse can be written as: where: is the pulse duration at the output of the dispersive medium and A’ the new amplitude. In the following, we shall refer to the Stokes and pump-and-probe pulses of typical two-color CARS implementations by using “Stokes” and “pump” subscripts in the relevant quantities. The pulses travel through several dispersing elements in the setup, including for example the scan and tube lenses, the microscope objective and any additional glass element introduced to tune their chirp. In the following, a will refer to the overall GDD. The amplitudes of the sum-frequency-generation signal A and of the four-wave mixing one A can be written as: where σ and σ describe the process efficiencies, and the “*” superscript indicates the complex conjugate of the field. Here positive time delays δt correspond to the pump beam arriving later than the Stokes beam, and we shall take into account only the contribution to the FWM signal generated from the absorption of two pump photons and the stimulated emission of a Stokes one. Second-harmonic generation (SHG) from the pump pulses and from the Stokes pulses will also be generated. Since these obviously do not depend on δt, they do not contain cross-correlation terms between the two different pulses and therefore will not be further considered in the following. The squared moduli of the Fourier transforms of A and A provide their spectra: The explicit expressions of the coefficients in Eq (5) are reported in S1 Table. From Eq (5) it is straightforward to notice that at fixed δt, the spectra are Gaussian curves with amplitudes (I and I) and centers (ω, and ω,) depending on δt. The amplitudes of the Gaussian curves as a function of δt are: where β and β are defined in S1 Table. The full-widths at half maximum (W) are: The FWHMs can be straightforwardly determined experimentally by monitoring I and I while moving the delay line. Inverting Eq (7) allows calculating the values of τout,pump and τout,Stokes: If the pulse durations are known, then the inversion of Eq (3) yields the GDDs at sample level. If the durations are not known but the center frequency of the SFG and FWM signals can be determined by acquiring their spectra, it is still possible to retrieve the values of the GDDs. To this end, it should be noted from Eq (5) that the center frequencies of the spectra shift linearly as a function of δt with rates (S): Inversion of Eq (9), after measuring the rates and calculating the values of τout,pump and τout,Stokes, yields the values of α and α: and, straightforwardly, of τ, τ, a, and a: yielding the complete characterization of the pulses at the sample level.

Materials and Methods

In the following, we shall show the validity of our approach by reporting the pulse characteristics of a CARS setup, schematically shown in Fig 1a, based on a fs laser for the generation of the pump-and-probe pulses with a wavelength of 810 nm and an OPO for the Stokes pulses (1060-nm wavelength). The beams propagate through two blocks of SF6 glass and two beam expanders. In the following we shall indicate with L (L) the length of the glass blocks placed in the pump (Stokes) path. Before the dichroic mirror D, the Stokes pulses are delayed by means of a motorized delay line. Photons reach and overfill the pupil of the microscope objective after being reflected by a pair of galvanometric mirrors and after propagating through a scan lens and a tube lens. Pump (Stokes) beam power measured after the objective is 20 mW (10 mW). Signals generated from the sample are collected by the microscope condenser, spectrally filtered to remove pump and Stokes photons, and analyzed by a spectrometer.
Fig 1

a) Schematic view of the CARS microscopy setup used in this work. Legend: fs laser (Chameleon Vision II). Is: Faraday optical isolator. OPO: Radiantis ORIA optical parametric oscillator. G, G’: dispersing SF6-glass blocks. BE: beam expanders. M: mirrors. DL: delay line. D: dichroic mirror. XY: galvanometric mirrors. Obj: Zeiss 32x N.A. = 0.85 C-Achroplan. S: sample. Cond: microscope condenser. BP: razor edge and band pass filters. Spectrometer: Jobin Yvon HR550 monochromator coupled to a Synapse CCD camera. The pump-and-probe (Stokes) beam path is shown as dotted (dashed) lines. b) Scanning-electron-microscope image of core/shell of the BaTiO3/Au nanoparticles used in this work (scale bar: 5 μm). Inset: magnified view of an individual nanoparticle showing decoration of the core with Au clusters (scale bar: 200 nm). c) Semi-log plot of the spectrum collected from a nanoparticle without G and G’ and with the pump-and-probe and Stokes beam overlapping in time.

a) Schematic view of the CARS microscopy setup used in this work. Legend: fs laser (Chameleon Vision II). Is: Faraday optical isolator. OPO: Radiantis ORIA optical parametric oscillator. G, G’: dispersing SF6-glass blocks. BE: beam expanders. M: mirrors. DL: delay line. D: dichroic mirror. XY: galvanometric mirrors. Obj: Zeiss 32x N.A. = 0.85 C-Achroplan. S: sample. Cond: microscope condenser. BP: razor edge and band pass filters. Spectrometer: Jobin Yvon HR550 monochromator coupled to a Synapse CCD camera. The pump-and-probe (Stokes) beam path is shown as dotted (dashed) lines. b) Scanning-electron-microscope image of core/shell of the BaTiO3/Au nanoparticles used in this work (scale bar: 5 μm). Inset: magnified view of an individual nanoparticle showing decoration of the core with Au clusters (scale bar: 200 nm). c) Semi-log plot of the spectrum collected from a nanoparticle without G and G’ and with the pump-and-probe and Stokes beam overlapping in time. As a sample material for the experiments, we chose core/shell BaTiO3/Au nanoparticles, (NPs), with an average diameter of 300 nm owing to their strong non-linear optical properties. A scanning-electron-microscope image of some NPs deposited on a microscope glass slide is shown in Fig 1b. The full synthesis and characterization details of the NPs are reported elsewhere [9,27].

Results and Discussion

In order to test our approach, we introduced known amounts of GDD in the pump and Stokes paths, determined the characteristics of the pulses, and checked them for consistency. The first test was carried out in the absence of SF6 glasses (L = L = 0) to evaluate the contribution of the set-up optics and sources alone to the total GDDs. The measured spectra as a function of the delay are reported in Fig 2a, where the top left (bottom left) panel shows the SFG (FWM) signal. Each spectrum was fitted with a Gaussian curve, yielding its center wavelength (shown in Fig 2a, left panels, as dashed curves, as a function of the delay) and its amplitude (shown in Fig 2b as open circles). The latter was fitted as a function of the delay again with Gaussian curves yielding W and W. Linear fits of the center wavelengths as a function of the delay in turns gives S and S. Eqs (8), (10) and (11) finally gives the pulse parameters, which were used to calculate the theoretical spectra (Fig 2a, right panels) for comparison with the experimental data.
Fig 2

(Color online) a) Measured (left panels) and simulated (right panels) sum-frequency generation (top row) and four-wave mixing (bottom row) signals with L = L = 0. Data and calculation results are shown normalized from 0 (dark) to 1 (white). The dashed lines represent the centers of the spectra determined by Gaussian fits of the data. b) Amplitude of the SFG (left) and FWM (right) spectra shown in Panel (a) as a function of the delay. Experimental data are shown as circles. Gaussian fits to the data are displayed as solid lines.

(Color online) a) Measured (left panels) and simulated (right panels) sum-frequency generation (top row) and four-wave mixing (bottom row) signals with L = L = 0. Data and calculation results are shown normalized from 0 (dark) to 1 (white). The dashed lines represent the centers of the spectra determined by Gaussian fits of the data. b) Amplitude of the SFG (left) and FWM (right) spectra shown in Panel (a) as a function of the delay. Experimental data are shown as circles. Gaussian fits to the data are displayed as solid lines. The parameters calculated for L = L = 0 and for the other configurations tested, averaged over 20 repetitions of each experiment, are reported in S2 Table and in S1 Fig. It is interesting to note that, without G and G’ (i. e. L = L = 0) even if both pulses cross similar amounts of glass in the setup, they reach the sample with a difference in GDD ~10000 fs2. This is quite likely due to a setting of the OPO cavity introducing a negative chirp in the Stokes beam that compensates the positive one introduced by the additional optics. We repeated the procedure after introducing L = 10 cm and L = 15 cm blocks of SF6 glass. The lengths were calculated in order to introduce approximately the same GDD for both the pump and Stokes beams, which is what one would be tempted to choose without knowing the real characteristics of the pulses. Calculating the refractive index of the SF6 glass by means of the Sellmeier equation (with parameters B1 = 1.72448482, B2 = 0.390104889, B3 = 1.04572858, C1 = 0.0134871947 μm2, C2 = 0.0569318095 μm2, C3 = 118.557185 μm2), the additional dispersions can be determined to be 19643 fs2 and 19746 fs2 for the pump and Stokes pulses respectively. The characterization of SFG and FWM, shown in Fig 3, yielded pulse durations–reported in S2 Table and in S1 Fig–consistent with the previous configurations and also GDD values in agreement with those determined for L = 0 and L = 0, increased by the calculated GDDs introduced by the glass. Similarly consistent parameters were obtained for L = 0 cm and L = 15 cm.
Fig 3

(Color online) a) Measured (left panels) and simulated (right panels) sum-frequency generation (top row) and four-wave mixing (bottom row) signals with L = 10 cm, L = 15 cm. Data and calculation results are shown normalized from 0 (dark) to 1 (white). The dashed lines represent the centers of the spectra determined by Gaussian fits of the data. b) Amplitude of the SFG (left) and FWM (right) spectra shown in Panel (a) as a function of the delay. Experimental data are shown as circles. Gaussian fits to the data are displayed as solid lines.

(Color online) a) Measured (left panels) and simulated (right panels) sum-frequency generation (top row) and four-wave mixing (bottom row) signals with L = 10 cm, L = 15 cm. Data and calculation results are shown normalized from 0 (dark) to 1 (white). The dashed lines represent the centers of the spectra determined by Gaussian fits of the data. b) Amplitude of the SFG (left) and FWM (right) spectra shown in Panel (a) as a function of the delay. Experimental data are shown as circles. Gaussian fits to the data are displayed as solid lines. The last combination of glass blocks used, namely L = 10 cm and L = 25 cm, was determined following the characterization in the other configurations in order to yield similar GDDs (~30000 fs2) for both pulses at the sample level. The increase in spectral resolution with the latter choice of glass lengths was demonstrated by measuring the spectrum of liquid methanol in the C-H stretch region (shown in Fig 4) and comparing it with the spectrum calculated using the pulse characteristics determined with our algorithm. The height of the Raman peaks and the amount of non-resonant background were adjusted to yield the best fit to the data, while their positions were fixed at 2944 cm-1 and 2836 cm-1, with widths of 34 and 21 cm-1 respectively, according to Ref. [28]. Fig 4 also shows the calculated spectra corresponding to L = 10 cm and L = 15 cm and to L = 0 and L = 0. The spectrum calculated with the pulse parameters determined with the analysis presented here provides a very good match to the experimental data, further confirming the validity of our approach.
Fig 4

(Color online) CARS spectra of liquid methanol.

The experimental spectrum, measured with L = 10 cm and L = 25 cm is shown with circles. The corresponding calculated spectrum is shown as a black solid line. Spectra calculated for L = 0 and L = 0 (dashed line) and L = 10 cm and L = 15 cm (dash-dotted line) are also reported for comparison.

(Color online) CARS spectra of liquid methanol.

The experimental spectrum, measured with L = 10 cm and L = 25 cm is shown with circles. The corresponding calculated spectrum is shown as a black solid line. Spectra calculated for L = 0 and L = 0 (dashed line) and L = 10 cm and L = 15 cm (dash-dotted line) are also reported for comparison.

Conclusions

We have shown that the pulse characteristics of a typical CARS setup can be easily retrieved, without using additional instrumentation, by analyzing the FWM and SFG signals from an appropriate non-linear sample. This information can be exploited to precisely determine and tune the chirp and pulse duration at sample level and therefore to optimize the performance of CARS microscopes in terms of spectral resolution and imaging contrast, and signal level for each specific application.

Summary of the pump-and-probe and Stokes pulse characteristics with different lengths of SF6 glass introduced in the beam paths.

The labels on the horizontal axis represent the lengths in cm of the glass introduced in the pump-and-probe path and in the Stokes one. Bars are pattern-filled depending on the amount of glass introduced in the pump-and-probe path (left column) or in the Stokes path (right column) to allow an easy comparison among the data. Error bars show the standard deviation calculated over 20 repetition of the acquisition. (EPS) Click here for additional data file.

Measured values of W as a function of the ratio between the Stokes and pump laser powers with L = 10 cm and L = 25 cm.

Within the measurement uncertainties, no variations of W were observed. In the case of probe-induced pump depletion, at large ratios a reduction of W is expected. (EPS) Click here for additional data file.

Testing of the beam profiles.

In order to verify the assumption of Gaussian beams, one of the beams was stretched with a block of glass and its temporal profile determined using the other pulse as a much narrower optical gate. The plot shows FWM intensities (open circles) as a function of δt with a 10-cm SF6 glass block in the Stokes-beam path. The red line shows the best Gaussian fit to the data (FWHM = 478 ± 3 fs) and the crosses displays the fit residuals. Crosses show the residuals, demonstrating the quality of the fit. Similarly, the pump beam was tested with analogous results. (EPS) Click here for additional data file.

Explicit expressions of the coefficients displayed in Eq (5).

(DOCX) Click here for additional data file.

Pulse parameters determined as a function of L and L.

Errors represent the standard deviation calculated over 20 repetition of the acquisition. (DOCX) Click here for additional data file.

Contribution to the spectral phase of the third-order dispersion term.

(DOCX) Click here for additional data file.
  15 in total

1.  Biophotonic probing of macromolecular transformations during apoptosis.

Authors:  Artem Pliss; Andrey N Kuzmin; Aliaksandr V Kachynski; Paras N Prasad
Journal:  Proc Natl Acad Sci U S A       Date:  2010-07-06       Impact factor: 11.205

2.  A large-field polarisation-resolved laser scanning microscope: applications to CARS imaging.

Authors:  G DE Vito; A Canta; P Marmiroli; V Piazza
Journal:  J Microsc       Date:  2015-07-30       Impact factor: 1.758

3.  Chemical imaging of tissue in vivo with video-rate coherent anti-Stokes Raman scattering microscopy.

Authors:  Conor L Evans; Eric O Potma; Mehron Puoris'haag; Daniel Côté; Charles P Lin; X Sunney Xie
Journal:  Proc Natl Acad Sci U S A       Date:  2005-11-01       Impact factor: 11.205

4.  Chirped coherent anti-stokes Raman scattering for high spectral resolution spectroscopy and chemically selective imaging.

Authors:  Kelly P Knutsen; Benjamin M Messer; Robert M Onorato; Richard J Saykally
Journal:  J Phys Chem B       Date:  2006-03-30       Impact factor: 2.991

5.  Single-shot measurement of the intensity and phase of an arbitrary ultrashort pulse by using frequency-resolved optical gating.

Authors:  D J Kane; R Trebino
Journal:  Opt Lett       Date:  1993-05-15       Impact factor: 3.776

6.  Cytocompatibility evaluation of gum Arabic-coated ultra-pure boron nitride nanotubes on human cells.

Authors:  Gianni Ciofani; Serena Del Turco; Antonella Rocca; Giuseppe de Vito; Valentina Cappello; Maho Yamaguchi; Xia Li; Barbara Mazzolai; Giuseppina Basta; Mauro Gemmi; Vincenzo Piazza; Dmitri Golberg; Virgilio Mattoli
Journal:  Nanomedicine (Lond)       Date:  2014-02-06       Impact factor: 5.307

7.  Rotating-polarization CARS microscopy: combining chemical and molecular orientation sensitivity.

Authors:  Giuseppe de Vito; Angelo Bifone; Vincenzo Piazza
Journal:  Opt Express       Date:  2012-12-31       Impact factor: 3.894

8.  Hyperspectral imaging with stimulated Raman scattering by chirped femtosecond lasers.

Authors:  Dan Fu; Gary Holtom; Christian Freudiger; Xu Zhang; Xiaoliang Sunney Xie
Journal:  J Phys Chem B       Date:  2013-01-18       Impact factor: 2.991

9.  Multimodal nonlinear optical imaging of collagen arrays.

Authors:  Christian P Pfeffer; Bjorn R Olsen; Feruz Ganikhanov; François Légaré
Journal:  J Struct Biol       Date:  2008-07-11       Impact factor: 2.867

10.  Barium titanate core--gold shell nanoparticles for hyperthermia treatments.

Authors:  Elmira FarrokhTakin; Gianni Ciofani; Gian Luigi Puleo; Giuseppe de Vito; Carlo Filippeschi; Barbara Mazzolai; Vincenzo Piazza; Virgilio Mattoli
Journal:  Int J Nanomedicine       Date:  2013-06-28
View more
  2 in total

1.  Lipid biochemical changes detected in normal appearing white matter of chronic multiple sclerosis by spectral coherent Raman imaging.

Authors:  K W C Poon; C Brideau; R Klaver; G J Schenk; J J Geurts; P K Stys
Journal:  Chem Sci       Date:  2018-01-02       Impact factor: 9.825

2.  Excitation parameters optimized for coherent anti-Stokes Raman scattering imaging of myelinated tissue.

Authors:  Craig Brideau; Kelvin W C Poon; Pina Colarusso; Peter K Stys
Journal:  J Biomed Opt       Date:  2019-04       Impact factor: 3.170

  2 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.